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(ABSTRACT) 

A study has been performed focusing on the calculation of sensitivities of dis- 
placements, velocities, accelerations, and stresses in linear, structural, transient 
response problems. One significant goal of the study was to develop and evaluate 
sensitivity calculation techniques suitable for large-order finite element analyses. 
Accordingly, approximation vectors such as vibration mode shapes are used to re- 
duce the dimensionality of the finite element model. Much of the research focused 
on the accuracy of both response quantities and sensitivities as a function of number 
of vectors used. 

Two types of sensitivity calculation techniques were developed and evaluated. 
The first type of technique is an overall finite difference method where the analysis 
is repeated for perturbed designs. The second type of technique is termed semi- 
analytical because it involves direct, analytical differentiation of the equations of 
motion with finite difference approximation of the coefficient matrices. To be com- 
putationally practical in large-order problems, the overall finite difference methods 
must use the approximation vectors from the original design in the analyses of the 
perturbed models. In several cases this fixed mode approach resulted in very poor 
approximations of the stress sensitivities. Almost all of the original modes were 
required for an accurate sensitivity and for small numbers of modes, the accuracy 
was extremely poor. To overcome this poor accuracy, two semi-analytical techniques 



were developed. The first technique accounts for the change in eigenvectors through 
approximate eigenvector derivatives. The second technique applies the mode accel- 
eration method of transient analysis to the sensitivity calculations. Both result in 
accurate values of the stress sensitivities with a small number of modes. In both 
techniques the computational cost is much less than would result if the vibration 
modes were recalculated and then used in an overall finite difference method. 
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Chapter 1 


Introduction 


1.1. Overview 

In the past ten years there has been increasing interest in calculating the deriva- 
tives of structural behavior with respect to problem parameters or design variables 
(i.e., sensitivities). One of the main uses of these sensitivities is in automated design 
procedures where a numerical algorithm is used to improve a structure by modify- 
ing the design parameters while satisfying prescribed constraints on the structural 
behavior. Most of the numerical algorithms used in these procedures require both 
an initial design and a set of sensitivities in order to decide how to improve the 
structure. Many references address this sensitivity calculation question within the 
context of automated structural design while others, such as this study, focus specifi- 
cally on issues related to the calculation of sensitivities. Other uses of sensitivities in 
structures problems include the system identification problem in structural dynam- 
ics and statistical structural analysis. References [1] and [2] provide a comprehensive 
review of work on calculating sensitivities in structural systems. 

It is clear from many references (e.g., ref. [1]) that most of the emphasis in 
structural optimization and the associated sensitivity calculation methods has been 
on static problems. This is not surprising since the majority of structural analyses 
themselves are static. The objective of the static analysis and sensitivity calculation 
problem, for linear systems, is to calculate the responses (e.g., displacements and 
stresses) and their derivatives with respect to structural parameters (e.g., member 
areas and thicknesses) which are assumed to be constant for all time. Techniques for 
both the analysis and sensitivity calculations have reached considerable maturity in 
the past ten to twenty years. 
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In many problems, however, the loading on the structure varies with time which 
causes the response of the structure also to vary as a function of time. Examples of 
such problems are a gust on an aircraft wing, an unbalanced engine in an automobile, 
or a building during an earthquake. In these cases it is important to predict stresses 
accurately as well as displacements (and possibly velocities and accelerations) as a 
function of time. Often it is sufficient to predict the maximum and minimum values 
of these response quantities. Similarly, the goal of the sensitivity analysis is the 
calculation of derivatives of these response quantities with respect to the structural 
parameters as a function of time or at the time points where the maximum or 
minimum responses occur. 

The introduction of the time parameter complicates the analysis in several 
ways. First, it changes the system of equations from a set of coupled algebraic 
equations to a set of coupled differential equations whose accurate solution may be 
difficult and computationally costly. Second, the amount of information that must 
be considered and evaluated to understand the response of the structure is increased 
by orders of magnitude. 

Most practiced static and dynamic analyses are currently performed using the 
finite element method. Since this technique replaces a continuum (infinite dimen- 
sional space) with a finite degree-of-freedom approximation, the question of required 
mesh refinement is a natural one. This is not an easy question to answer because the 
convergence of the approximation as the mesh is refined depends on the quantity 
being considered. Usually, the fundamental unknowns are the displacements and 
rotations at the finite element nodes. In theses cases, the convergence of derivatives 
of displacements with respect to a spatial parameter (stresses), with respect to time 
(velocities, accelerations), or with respect to a structural parameter (sensitivities) 
will be worse than the convergence of the displacements themselves. 
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After the structure has been discretized using the finite element method, yet 
another approximation is usually introduced in linear dynamics problems. The be- 
havior of the structure is represented by a reduced set of basis functions (frequently 
natural vibration modes) in order to simplify the solution of the transient response 
problem. This approximation introduces another set of concerns over accuracy of 
the response quantities and their sensitivities. 

Other errors in transient analysis or sensitivity calculations, which rarely occur 
in static analyses, are due to the truncation error of finite difference operators. 
This problem occurs with the use of numerical integration techniques in solving the 
coupled differential equations in the transient problem. This problem also occurs 
when difference approximations are used in the calculation of sensitivities. Roundoff 
errors, due to the finite precision arithmetic on digital computers, are also more of 
a concern in transient or sensitivity analyses than in simple, static analyses. 

All of the above discussed complexities in transient analysis coupled with the 
problems of sensitivity analysis have slowed the progress in the development of sensi- 
tivity calculation techniques for transient response problems. However, substantial 
progress has been made. Some of the important, previous work in optimization of 
structures under transient loads and calculation of sensitivities in transient response 
problems is discussed below. 

1.2. Review of Previous Pertinent Work 

Reference [3] is one of the earliest papers dealing with optimization of struc- 
tures under transient loads. In this work Fox and Kapoor consider the minimum 
mass design of frame structures under an applied base motion subject to constraints 
on deflections and stresses. The equations of motion are uncoupled using vibration 
modes and solved for the maximum value of the modal response using a shock 
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spectrum approach. A considerable simplification is introduced by directly sum- 
ming the maximum modal responses, thereby removing time as a parameter in the 
calculations. 

In references [4] and [5] Cassis and Schmit present procedures for the automated 
design of plane frames under general transient loading. The dynamic analysis is 
performed using modal superposition and only modal damping is allowed. Integral 
forms of the time dependent constraints are used. Sensitivities are calculated using 
an explicit differentiation of the dynamic equations along with exact calculation 
of the required eigenvalue and eigenvector derivatives. Effects of finite element 
discretization and modal truncation on the sensitivities or final, optimized designs 
were not considered. 

In the past 10 years, other researchers have considered the application of gen- 
eral sensitivity theory to the problem of dynamic mechanical systems. Reference 
[1] summarizes this work and describes three basic approaches which have been 
employed. In the first method, called the direct method, the equations of motion 
are directly differentiated and solved. A second method offers the advantage of 
reduced computational cost when there are more design variables than constraints 
on response quantities. In this method, called the adjoint method, the sensitivity 
equations are rewritten in terms of a newly defined adjoint vector. After solving 
this new system for the adjoint vector, the calculation of the sensitivities of the re- 
sponse constraints with respect to each of the design variables is straightforward. In 
the third method, called the Green’s function method, the derivatives are obtained 
in terms of the Green’s function of the equations of motion. Although the results 
from all three methods are theoretically identical, their relative computational effi- 
ciency depends on the relative numbers of design variables, degrees of freedom, and 
constraints. 
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Haug, Arora and their co-workers have made considerable progress in address- 
ing many of the problems in the optimal design of mechanical systems under dy- 
namic loadings. Much of their early work was spent studying a “state-space” or 
adjoint variable approach to calculating sensitivities. References [6], [7], [8], and [9] 
should be noted. These references consider application to both elastic structured 
design and machine design problems that often have the additional complexity of 
nonlinear equations of motion. However, most of these examples have involved few 
degrees-of-freedom or design variables. A more recent paper by Haug ([10]) ex- 
tended the sensitivity analyses of previous papers to include additional, algebraic 
constraint equations that are often present in machine design problems. Also, sen- 
sitivity equations for second derivatives are presented. 

The adjoint method is particularly attractive when a transient constraint is 
integrated over time to produce a single constraint because the total number of 
constraints is often small. However, the loss of information in this integral formu- 
lation and its disadvantages are noted in reference [11]. Given the danger of having 
only a single “worst-case” value of the constraint function in time, reference [11] 
proposed including all local maximum points of the constraint function in the con- 
straint set. A significant disadvantage of this approach is that for “jagged” response 
functions, there can be a large number of redundant, local maxima. This important 
problem of constraint definition was also considered in reference [8], where several 
methods for obtaining a few, important constraints at discrete points in time were 
proposed. 

Both direct and adjoint sensitivity methods for a nonlinear, hysteretic structure 
are presented in reference [12]. Because of the nonlinearities, numerical integration 
of the full coupled system is required. 

A recent approach in sensitivity analysis has been to write sensitivity expres- 
sions for the solid continuum prior to discretizing the system. This approach is 


5 



especially attractive when shape-type design variables are being considered because 
the design variable itself often represents a continuous region on the surface of the 
body. Reference [13] uses the concept of the material derivative to calculate shape 
derivatives of a continuum under dynamic loads. In reference [14] expressions for 
shape sensitivities of a continuum considering material nonlinearities and dynamic 
effects are written using a variational approach. 

1.3. Objectives and Scope 

The purpose of the study reported herein is to investigate methods for calcu- 
lating sensitivities in linear transient structural response problems. Very general 
forms of external loading on the structure and damping are permitted. In any nu- 
merical algorithm, both accuracy and computational efficiency are concerns. Errors 
in the sensitivities due to factors such as the finite element mesh, truncation of the 
basis vector set in the transient analysis, and finite difference approximations in 
the sensitivity and numerical integration procedures are considered. An objective 
of the study is to identify approaches to sensitivity analysis that are appropriate 
for large-scale structural analysis. This is emphasized in the selection of the algo- 
rithms and in a study of the relative computational efficiency of several competing 
methods. 

Three transient response problems are considered in detail- a five-span, simply 
supported beam, a composite aircraft wing, and a cantilever beam with a cross 
section that varies along its length. None of these three problems is large. However, 
each problem includes ingredients which make the sensitivity analysis computation- 
ally difficult. 
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Chapter 2 


Equations of Motion and Solution 
2.1. Governing Equations 

The equations of motion for a damped, linear structural system can be written 
as 

Mu + Cu + Ku = p(f) (2-1) 

which is a set of n 9 , coupled differential equations. M, C, and K are the system 
mass, damping, and stiffness matrices, respectively. Frequently it is possible to 
separate the loading vector, p, into a product of a vector describing the spatial 
distribution of the loading, f, and a scalar function of time g(t) as 

P (i) = 9{t)f (2-2) 

Often equations (2.1) are the result of a large finite element model and are 
therefore of large order. One way to characterize the behavior of this system is by 
examining the eigenvalues of the undamped system 


K<f> — u>jM^ = 0, j=l,...,n, (2.3) 

For most large, structural systems Eqs. (2.1) are “stiff”; the condition number, 
is many orders of magnitude. 

The external loading also has a major effect on the dynamic response of the sys- 
tem. Impulsive loads where g(t) changes rapidly relative to the periods associated 
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with the smallest u>j tend to produce a response history with significant high fre- 
quency components. Loads that are applied slowly relative to the vibration periods 
of the u ij produce a predominantly low frequency response history. 

Two basic approaches are available for the solution of equations (2.1). The 
first is to numerically integrate the equations in a step-by-step manner. In implicit 
integration techniques, the time step must be a fraction of the period associated 
with the largest u>j significantly excited by the loading in order to obtain an accurate 
solution. The well-known Newmark method (see for example ref. [15]) is an example 
of such an integration technique. In explicit integration techniques the time step 
must be a fraction of the period associated with u ng in order for the solution process 
to be numerically stable. Using either technique, the computational work is large 
because equations (2.1) are of large order. 

An alternative to directly solving equations (2.1) is to solve an approximate, 
reduced order problem instead. This is the preferred approach for most linear 
structural dynamics problems. The details of the techniques used to reduce the 
order of the dynamic system are discussed below. 

2.2. Reduction Techniques 

The first step in applying a reduction technique to the solution of equations 
(2.1) is to approximate the solution by n r basis functions 


u = $q (2-4) 

where n r is usually much less than n g . Then a reduced set of equations can be 
written 


Mq+Cq + Rq = g{t)l 


(2.5) 
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where 


M = * t M# 


e = $ T c$ 

R = i r K$ 

r = * T f 


( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 
(2.9) 


If the number of vectors in is equal to the size of the original system, n g , and the 
vectors in 4* are linearly independent, the transformation of equation (2.4) is exact. 
Usually, though, n r « n g and the solution to the full system (equations (2.1)) is 
only approximated by the solution to the reduced system (equations (2.5)). The 
quality of this approximation as the number of vectors in # is increased is a key 
concern in evaluating the effectiveness of a particular reduction technique. 

In all of the reduction methods considered herein, the first n r vectors of the 
set are taken as the reduced basis. Alternate approaches are available for assessing 
the importance of a given vector prior to solution of the reduced system and then 
discarding the vector if its contribution is insignificant. These approaches are not 
considered here because the cost of generating the set of vectors $ is often high and 
the cost of solving equations 2.5 is often fairly low. 

2.2.1. Mode Displacement Method 

The most widely used reduction technique is the traditional mode displacement 
method. In this method, equations (2.3) are solved for the set of vibration modes 
with lowest n r frequencies and modes. This set of vibration modes is used as the set 
of basis functions *&. When the system is undamped (C = 0 in equations (2.1)) or 
C can be expressed as a linear combination of M and K, equations (2.5) represent 
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a set of uncoupled differential equations which can be solved independently. If the 
eigenvectors are scaled so that <f>J = 1, the uncoupled equations can be written 

q% + 2 &u>i 9 i + u>?qi = g{t)fi, i = l,...,n r (2.10) 

where & is the modal damping ratio. For certain forms of external loading, such as 
g(t) represented as a piecewise linear function of time, an exact, explicit solution is 
available. This approach is described in [16] and is used in NASTRAN [17]. 

Equations (2.1) are the result of a given finite element approximation designed 
to model the behavior of the dynamic system. The goal of the reduction methods 
discussed in this section is to achieve an accurate approximation to the solution 
to equations (2.1) using a small number of basis vectors. As discussed above, the 
vibration modes are the most commonly used basis functions in linear structured 
dynamics. There are two cases, however, where a large number of modes are re- 
quired for an accurate solution of equations (2.1) and therefore the performance of 
the mode displacement method is poor. In the first case, if the structure is loaded 
in an impulsive manner, many high frequency modes tend to be excited. These 
high frequency modes must be included in the analysis since their contribution to 
the total response is significant. In the second case, if the response of the structure 
contains a large static component, the linear combination of vibration modes can 
do a poor job of approximating the static deflection shape. The reduction methods 
discussed below alleviate this second accuracy problem with the mode displacement 
method. 
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2.2.2. Mode Acceleration Method 


To alleviate the poor accuracy of the mode displacement method due to its poor 
representation of the static component in the response, a method was proposed by 
Williams (ref. [18]) called the mode acceleration method which is described in its 
modern computational forms in refs. [16] and [19]. The mode acceleration method 
can be derived by rewriting equations (2.1) as 


u(<) = g(t)K- 1 f-K- 1 Cu-K" 1 Mu (2.11) 


The first term in equations (2.11) is the quasi-static solution with load amplitude 
determined by g(t). This term is calculated by solving the equations 

Ku, = f (2.12) 


This solution is carried out in the standard way by first factoring K into a product 
of upper and lower triangular matrices and then performing a forward and backward 
substitution operation to obtain u t . The other two terms are calculated using the 
solution u and u from the mode displacement solution. In these terms K 1 is 
calculated as follows. Equations (2.3) can be rewritten in matrix form as 

K$n -2 = M$ (2.13) 


where here # is the full set of n g eigenvectors and fl 2 is 


rl/u,? 



With the eigenvectors scaled so that 


= I 


(2.14) 


(2.15) 


11 




equation (2.13) can be written as 


$ T K*n- 2 = I (2.16) 

Premultiplying by (3» T ) -1 and postmultiplying by yields 

K$0- 2 $ t = I (2.17) 

or 

K 1 = $n- 2 $ T (2.18) 

When $ contains less than the full n g eigenvectors, this expression for K" 1 is only 
approximate. However, since ii is obtained from the mode displacement solution 
based on n r modes ($q), K _1 Mu is exactly equal to $0 -2 q and no approxima- 
tion results from introducing equations (2.18). For the damping term, introducing 
equations (2.18) with n r vectors in $ is not exact. However, this is a convenient 
approximation especially when modal damping is used. Consistent with these con- 
siderations, equations (2.11) can be rewritten as 

u(t ) = flr(i)K -1 f - $n" 2 Cq - $n~ 2 q (2.19) 

The key to the effectiveness of this method is that the static solution is included 

explicitly in the solution. It is also simple to apply since it essentially just super- 
imposes the static and mode displacement solutions. Since q and q are obtained 
from the mode displacement solution, u and u are identical to the values obtained 
in the mode displacement method. 
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2.2.3. Static Mode Method 


An alternative approach to the mode acceleration method that accounts for 
the static solution slightly differently is termed the static mode method herein. In 
this method, the static solution is included as an additional “mode” in forming 
the reduced equations (2.5). The procedure begins by calculating a set of n r — 1 
eigenvectors ^ using equations (2.3). Then the static solution is calculated as 

K<^ = f (2.20) 

To improve the orthogonality of the basis vectors, the components of the vibration 
modes are removed from the static solution using the Gram-Schmidt process 

4>i = — $c (2.21) 

where 

c = (2.22) 

The vector <f>\ is then concatenated with $ to yield a new which is the complete 
basis. Equations (2.5) now become coupled and can be solved directly or reduced 
to an uncoupled form using the following procedure. First the reduced eigenvalue 
problem 


MZA + RZ = 0 (2.23) 

is solved for the n r x n r diagonal matrix of eigenvalues, A, and the n r X n r matrix 
of eigenvectors Z. Now a new set of basis vectors can be written 


$ = $Z 


(2.24) 
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When # is substituted for $ in equations, (2.6), (2.7), (2.8), and (2.9), an uncoupled 
system results when C is of the special form described in section 2.2.1. 

The static mode method is similar to the mode acceleration method in that 
the static displacement vector is explicitly included in the solution. However, in 
the mode acceleration method the amplitude of the static displacement vector is 
not an unknown but is determined by g(t) while in the static mode method the 
amplitude varies to possibly improve the solution. Also, this static displacement 
vector participates in the calculation of u and ii to possibly improve them as well. 

2.2.4. Ritz-Wilson-Lanczos Method 


A fourth method, which has become popular in the past few years is termed the 
Ritz-Wilson-Lanczos (RWL) method and is described in refs. [20], [21], and [22]. 
Instead of using eigenvectors of the structure, this method uses a set of Lanczos 
vectors to form the reduced equations. The algorithm used here follows that in 
reference [20]. The first vector is obtained by solving the static equations (2.20) 
and then scaling so that 


4>i = (2-25) 

The vectors i — 2, . . . ,m are obtained as follows. First, 


K<^ = (2.26) 

is solved for Then 4 > i is made M-orthogonal with respect to all previously 
generated vectors using a Gram-Schmidt process 
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(2.27) 


»-i 

4>i = 4>i -^2 Ci ^j 

3 = 1 

where 


°ij - 4 >J M 4>i 


(2.28) 


scaling gives 


4>i = 4>i/(4>J (2.29) 

It has been pointed out in many references (e.g., [21]) that the M-orthogonalization 
(eq. (2.27)) is theoretically required only with respect to the two previously com- 
puted vectors. However, it is also well known that roundoff errors cause the Lanczos 
vectors to become less and less orthogonal. Performing the Gram-Schmidt operation 
with respect to all previously generated vectors will not insure the M-orthogonality 
of the vectors. However, it can improve the orthogonality in some cases. 

Following reference [20], a final step is performed in generating the basis vectors 
to produce an uncoupled dynamic system. Of course this is useful only when the 
system is undamped or C is assumed to be diagonal. As in the static mode method 
a reduced-order eigenvalue problem is solved (eq. (2.23)) and a new set of basis 
vectors produced. The process of explicitly computing the reduced stiffness and 
mass matrices required in equations (2.23) helps alleviate the problems caused by 
the lack of orthogonality of the Lanczos vectors. The matrices M and R in (2.23) 
are assumed to be full. That is, no assumptions are made that particular terms in 
M and R are zero based on the properties of the vectors. 
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2.3. Transient Response Solution Method 


When the above reduction methods are used and general damping is included 
in the model, equations (2.5) are coupled. In principle any of the implicit or explicit 
numerical integration methods used for solving equations (2.1) could be used to solve 
equations (2.5). In contrast to equations (2.1), however, equations (2.5) are low 
order, not stiff, and the primary concern is accurately integrating every equation in 
the system. Therefore an integration method which reduces truncation errors in the 
solution is highly desirable. Accuracy is especially important in sensitivity analyses 
because errors in the solution process are usually magnified in the calculation of 
derivatives. 

An approach that allows the use of moderately large time steps and makes the 
truncation error very small is called the matrix series expansion method in [23] and 
the transfer matrix method in [24] and [25] when applied to structural dynamics 
problems, and is often referred to as a Taylor series method in numerical analysis 
texts (e.g., ref. [26]). This method expands the solution in a Taylor series where the 
number of terms determines the accuracy of the approximation. Using this series 
an expression can be written for the solution at time t + At in terms of the solution 
and load at time t. 


Jq(t + A<)\ _[W„ W 12 l/q(t)\ [N n N 12 f g(t)f 1 
\q(f + A<)J [W 21 W 22 J \q(t)J + [N 21 N 22 \ <7(f)f J 


(2.30) 


It has been assumed here that the time variation of the load g(t) is approximated 
as a piecewise linear function of time and therefore the second and higher order 
derivatives equal zero. Expressions for W,y and N tJ are fairly complex and can be 
found in references [23] and [24]. The values of the coefficients W ij and Njj depend 
on the number of terms taken in the series. 
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The convergence properties of the Wjj series for an undamped, single degree- 
of- freedom case can be studied by considering the following Taylor series expansion 


cos u At — 1 — 


(u>At) 2 


(a>A<) 4 

4! 


(u>At) 6 (u>A <) 8 
~ 6 ! + 8 ! 


(2.31) 


It is well known that roundoff errors due to finite precision arithmetic will cause 
large errors in this series for “large” values of uAt. Thus if u; is taken as u? nr , 
an upper bound on At can be estimated based on roundoff error. In practice At 
will usually need to be much smaller than this upper bound value for two reasons. 
The first reason is that the input load history may be a complicated function of 
time and At must be small enough to accurately sample this loading. The second, 
more important reason, is that the At must be small enough to accurately sample 
the history of the output quantities. If At is larger than the smallest significant 
period of response, peak values of the response quantities will likely be missed. 
Accordingly, in the studies reported herein At was taken to be approximately 1/8 
of the smallest period. Since the number of terms in the series has only a very small 
effect on the computational cost of the method, 50 terms were used in this study 
to make the truncations errors negligible. 
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Chapter 3 


Critical Point Constraint 


3.1. Constraint Formulation 

The general form of the constraint equation is 


9i(x,t) < 0 (3.1) 

where x is a vector of design variables and i is time. An effective approach for 
insuring that this constraint is satisfied for all values of t is the “critical point 
constraint” approach described in ref. [27] pages 168-169. In this approach a set of 
peak values of the function gi (denoted critical points) is selected. An obvious point 
to include is the time with the “worst” value of g However, if only this point is 
included, an optimization process modifying a structure based on this information 
might unknowingly produce a design where the constraint is violated at another time 
point. To guard against this possibility, a number of important peaks are selected. 
References [28] and [29] consider in detail the efficient location of critical points 
in large-scale structures problems with many constraints. This chapter presents a 
method for selecting the most important peaks as critical points. 

In the work reported herein, constraints are assumed to be placed on the dis- 
placements, velocities, accelerations, and stresses in the structure. All of these 
constraints are treated similarly. Thus the critical point constraint formulation will 
be illustrated for the case of displacements. Constraints are placed on selected 
displacements such that 

|«i(x,f)| ^ U a l low (3.2) 
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where Ui axe the displacements at specific points in the structure and u a u ow is 
the absolute maximum allowable value of the displacement. The critical values of 
this constraint occur at points in time where U{ has the largest magnitude. These 
are identified by examining every value of u along the response history. In the 
implementation here, each constraint is assumed to have a specified number of 
critical points; five critical points for each Uj are selected. Values of u where du/dt = 
0 or values of u at the end points of the time interval are local maxima of gi and 
are termed candidate critical points. 

3.2. Selection of Critical Points 

The procedure for selecting the critical points from these candidates can best be 
explained by referring to an example displacement time history shown in figure 3.1. 
The critical points are labeled with circled numbers and a few of the many candidate 
critical points are labeled with circled letters. The selection criteria applied to every 
candidate critical point will be explained by considering these few candidate points. 
Candidate critical points a and c were discarded because the absolute values of 
the displacements at these points were smaller than those at the five other critical 
points. The criterion for discarding candidate points b, d, and e is slightly more 
complicated. From figure 3.1 it can be seen that all three candidate points have 
larger displacement magnitudes than that of critical point 1, for example. However, 
candidate points b, d, and e are all part of “major” peaks where a critical point 
is selected. A second criterion applied to the selection process is a requirement 
that only one critical point from each major peak be selected. This insures that 
the critical points represent the total dynamic response rather than just the high 
frequency undulations on, at worst, a single major peak. 

A major peak is identified with the following procedure. Whenever a critical 
point is selected after comparing its magnitude with that at other critical points, a 
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Figure 3.1. Example displacement time history illustrating the critical point con- 
straint selection process. 

special screening process is activated. This screening process tests the displacement 
at every subsequent time point to determine if it differs from that at this last selected 
critical point by at least a specified percentage (25% for the studies reported herein ). 
If so, all subsequent time points are no longer considered part of the current major 
peak. Any candidate critical points identified while this special screening process is 
in effect are compared only against the last selected critical point. 

An example is the major peak in figure 3.1 which contains points d and 4. In 
the selection process, point d is initially selected as a critical point and the screening 
process is activated. The three points where du/dt = 0 between point d and point 
4 are recognized to be part of the same major peak as d, but since the magnitude 
of the displacements at these points is smaller than at point d, they are discarded. 
Point 4 is also part of the same major peak as point d but since the displacement 
magnitude there is larger than at point d, it replaces point d as a critical point. 
Before the next candidate critical point is considered, the displacement has changed 
from that at point 4 by more than 25% and so is considered to be on a new major 
peak. 
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3.3. Derivatives of Critical Point Constraints 


Once the critical points have been identified for the nominal design, these can 
used in calculating sensitivities. Reference [27] demonstrates that the change in 
time location of critical points can be neglected in calculating derivatives of peak 
values with respect to design variables by examining the expression for the total 
derivative of gi with respect to a design variable x. Considering a constraint g(x,t) 
at a critical time t c 


dg(x,t c ) dg dg dt c 

dx ~0x dt dx ’ 

The last term in equation (3.3) is always zero because at interior critical points 
dg/dt = 0 and on the boundary dt c /dx = 0. Accordingly, the sensitivity calcu- 
lations need to be performed only at the specific times where the critical points 
have been identified. This can result in a considerable savings in computational 
time, especially when there are many constraints, many time points, or many basis 
functions used to represent the response. The details of each sensitivity calculation 
method are discussed in the next chapter. 
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Chapter 4 


Methods For Calculating Sensitivities 
4.1. Finite Difference Methods 

Both the forward difference and central difference methods have been used in 
this study to calculate sensitivities. The well known forward difference approxima- 
tion to du/dx , 


A u u(x + Ax) — u(x) 

Ax Ax 

and central difference approximation 


(4.1) 


Au u(x + Ax) — u(x — Ax) 
Ax 2Ax 


(4.2) 


are used. The truncation error for the forward difference approximation is 


ex{ Ax) 


Ax d 2 u 

~Ylx^ 


(x + £Ax) 


o < c < i 


(4.3) 


and is 


e r (Ax) = + C Ax ) 0<C<1 (4.4) 

for the central difference approximation. 

In applying equations (4.1) and (4.2), the selection of difference step size Ax 
is a concern. Selection of a large step size results in errors in the derivative due 
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to truncation of the operator (eqs. (4.3) and (4.4)). Selection of a small step size 
can lead to errors in the derivative due to the limited floating point precision of 
the computer or algorithmic inaccuracies in calculating u (condition errors). It is 
not uncommon with the forward difference method (eq. (4.1)) that no acceptable 
value exists for Aa; to produce an accurate value of du/dx considering the conflicting 
requirements of minimizing truncation and condition errors. Because the truncation 
error associated with the operator of equation (4.2) is typically less than that of 
equation (4.1), it is possible to use a larger finite difference step size. The larger 
Ax reduces the condition error from the function evaluations and results in a more 
accurate value of du/dx. However, the necessity of two function evaluations needed 
for equation (4.2) makes the procedure computationally more costly. 

4.1.1. Using Vibration Modes as Basis Functions 

For many of the studies herein the natural vibration mode shapes are used as 
basis functions to represent the transient response. In calculating the response of 
the perturbed design in equation (4.1) and the two perturbed designs in equation 
(4.2) some computational savings are possible relative to the computations for the 
initial design. 

If the mode shapes for the initial design are used to represent the perturbed 
design, the cost of re-solving the eigenvalue problem is eliminated. However, the 
reduced set of equations for the perturbed system must still be formed and M, C, K 
are now full. This coupled system is then solved using the matrix series expansion 
method described in section 2.3. 

If the updated mode shapes for the perturbed design are used in the analysis, 
many eigensolution procedures, such as the subspace iteration used here, can begin 
with the mode shapes from the initial design as approximations. Since the pertur- 
bation in the design is small, the subspace iteration procedure converges rapidly. 
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However, at least one factorization of K is required. For large finite element models 
this can be the largest part of the computational cost. For most of the studies in 
this paper, the forward difference method used the initial mode shapes to represent 
the perturbed design. Because the central difference method was used for reference 
values of derivatives, updated mode shapes were calculated for the two required per- 
turbed designs. In both cases, because of the critical point constraint formulation, 
the transformation from modal coordinates to physical coordinates (displacements, 
stresses, etc.) is performed only at the critical times instead of at all time points. 

4.2. Semi-Analytical Methods 

The direct method for sensitivity calculation is derived by differentiating equa- 
tion (2.1). The derivation presented here follows that on pages 169-171 in ref. [27]. 
After differentiating equations (2.1) with respect to the design variable, x , the result 
is 


m , du ^ du du df . , dM „ dC . dK 

dx dx dx dx yK ’ dx dx dx 


(4.5) 


This system of differential equations of order n g could be solved directly for the 
sensitivities, du/dx, dii/dx, and du/dx. However, just as for the response equations, 
it is more efficient to consider a reduced form of the sensitivity equations which can 
be obtained by differentiating equations (2.5) with respect to x yielding 

dR 


- dq _ dq - dq d? dM 

M-2 + CJ + K-3 = — 9 (() _ —q 
ax ax ax ax ax 


dC . 
d^ q 


dx 


(4.6) 


The first step in forming this equation is the calculation of the derivatives of f, M, 
C, and R (equations (2.9), (2. 6), (2. 7), (2.8)) with respect to x. Using equations (2.9) 
the derivative of ? with respect to x can be written as 


df d$ T , ~df 
— = f + $ — 

dx dx dx 


(4.7) 
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The force f is frequently not a function of the design variables which simplifies 
equation (4.7). Also using equation (2.6) the derivative of M with respect to z can 
be written as 


dM 

dx 


$ T - — $ H — M# + # t M- 


dx 


dx 


dx 


(4.8) 


Similar expressions can be written for the derivatives of C and R. 

The derivative dM/dx in equation (4.8) (similarly for dC/dx and dK/dz) is 
in general difficult to calculate because the finite element model may be composed 
of diverse element types whose properties are complicated functions of the design 
variables x. For this reason, these derivatives are often replaced with finite differ- 
ence approximations. This combination of analytical differentiation of the response 
equations with finite difference matrix derivatives is known as a semi-analytical 
approach. The semi- analytical methods presented herein for calculating transient 
response quantities all use the forward difference operator to approximate dM / dx , 
dC/dz, dK/dz, and df/dz. For several important classes of design variables, how- 
ever, M, C, and K are linear functions of z. For example, M and K in a finite 
element model composed of truss members are linear functions of member cross 
sectional area. In these cases there are no truncation errors and large finite dif- 
ference step sizes can be used to reduce the condition error and produce accurate 
derivatives. 


Calculation of the first term in equation (4.7) and the second and third terms in 
equation (4.8) depends on the particular choice of basis functions $. Considerable 
reduction in computational cost results if the vectors $ are taken to be independent 
of z, that is fixed. Methods are also available to approximate d4?/dz which are less 
costly than exact methods. Two semi- analytical procedures which address these 
concerns are discussed below. 


25 


4.2.1. Fixed Mode Semi- Analytical Formulation 


If the basis vectors are assumed not be functions of the design variables z, 
d&/dx = 0. This significantly simplifies equations (4.7) and (4.8). After forming 
the derivatives of ?, M, C, and K, the right-hand side of equations (4.6) can be 
formed using q , q, and q from the solution of equations (2.5). The matrix series 
expansion method insures that accurate values of q, q, and q are available for this 
step. Equations (4.6) can then be integrated to yield dq/dz, dq/dz, and dq/dz. 
This fixed mode, semi-analytical implementation of the direct method will be called 
just the semi- analytical method herein. 

4.2.2. Variable Mode Semi-Analytical Formulation 

If the basis functions are assumed to be functions of z, the calculation of 
d#/dz either exactly or approximately is required to form equations (4.7) and (4.8). 
Vibration modes are the most popular basis functions and the calculation of their 
derivatives has been studied extensively. Reference [30], for example, surveys several 
methods for calculating derivatives of vibration modes shapes from a computational 
point of view. One of the most popular methods, Nelson’s method (ref. [31]), 
requires a factorization of the system equations for each eigenvector considered. 
This can be a considerable computational burden for large systems. Since the 
overall objective here is the accurate calculation of transient response sensitivities, 
not eigenvector sensitivities, it seems desirable to investigate cheaper, approximate 
methods for calculating d$/dz. 

One approximate method for calculating the eigenvector derivatives is similar 
to the modal approach for transient analysis. This modal method approximates 
each eigenvector derivative as a linear combination of the modes themselves. In 
many cases, however, a very large number of eigenvectors are required for accurate 
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derivatives. Furthermore, the eigenvector derivative approximation produced by 
this method can’t improve the transient response sensitivities because they are 
contained entirely within the span of the modes themselves. 


A method proposed by Wang (ref. [32]) to alleviate the poor performance of 
the modal method also improves the transient sensitivities. This modified modal 
method is derived by first differentiating equations (2.3) to yield 


(4.9) 


This equation can not be solved directly since the left-hand-side is singular. Wang’s 
approach, however, was to calculate a pseudo-static solution to this equation by 
neglecting u>*M on the left-hand-side of equations (4.9). The solution to this pseudo- 
static equation introduces the change in basis associated with changes in the design 
variables and is significant in improving the transient response sensitivities. The 
mode shape derivative can then be written as 


d*j 

dx 


(£).*!*••• 


(4.10) 


where ( d$j/dx) a is the pseudo-static contribution. The coefficients, Ajk are ob- 
tained by substituting equation (4.10) into equation (4.9), multiplying by and 
simplifying as 


Ajk 


/j 2 7' ( dK. 2 dM \ ^ . 

1 dx dx ) 9 3 


- w fc) 


k±j 


or 


dM 




k = j 


(4.11) 


(4.12) 


Given these approximate values of eigenvector derivatives, equations (4.7) and (4.8) 
can be formed. Then equation (4.6) can be solved for dq/dx, dq/dx, and dq/dx 
just as in the fixed mode semi- analytical method. 
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4.2.3. Recovery of Physical Sensitivities 


Given dq/dx, the derivative of the physical displacement vector du/dx can be 
written as 

Att Art AS> 

(4.13) 


du _ dq d$ 

dx dx dx H 


with similar expressions for du/dx and du/dx. The calculation of stresses begins 
with 

cr = Su (4.14) 

where S is the stress transformation matrix. Substituting eq. (2.4) yields 


cr = S$q 


(4.15) 


Differentiating eq. (4.15), the stress sensitivities can be written as 


dcr 

dx 


_ ^ dq dS _ _ d$ 

ax uX dx 


(4.16) 


The matrix dS/dx is approximated using the forward difference operator. Because 
of the critical point constraint formulation, the transformation from these modal 


quantities to physical displacements, velocities, accelerations, and stresses is per- 


formed only at the critical times. 


4.2.4. Mode Acceleration Method 


The mode acceleration method was presented in Chapter 2 as a technique for 
improving the dynamic displacements and stresses when the static component is 
significant. It is also possible that it can improve the sensitivities of displacements 
and stresses. An expression for the sensitivities using the mode acceleration method 
is obtained by first rearranging equations (4.5) to yield 
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(4.17) 


du(t) _ j f df . . dK _ dii dC . _ , du dM _ 

dx [ dx dx dx dx dx dx 

If a reduced basis approximation is applied uniformly to every term in equations 
(4.17), the resulting du/dx would agree with that obtained from the solution of 
equation (4.6). The objective of a mode acceleration solution is to selectively apply 
the modal approximation to equation (4.17) with the goal of improving the values 
of du/dx. In applying the mode acceleration method to the transient response 
problem (eq. (2.11)) u and u are obtained from the mode displacement method. 
Here, in applying the mode acceleration method to the sensitivity equations, u is 
obtained from the solution to equations (2.19) and du/dx and du/dx are obtained 
from the solution to the mode displacement, semi-analytical equation ((4.6)). In 
the derivation here, the modes $ are assumed to be fixed. Substituting equations 
(2.19) and u = «fcq, u = 3>q, du/dx — #dq/dx, and du/dx = &dq/dx into (4.17) 
yields 



The modal approximation for K -1 (equation (2.18)) is introduced into all terms 
in equation (4.18) that involve damping just as in the mode acceleration solution 
described in Chapter 2. It was also pointed out in Chapter 2 that K -1 Mi£ in 
equation (4.18) is exactly $n -2 . Based on these considerations, equation (4.18) 
can be simplified to yield the mode acceleration solution of the sensitivity equations 
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du(t ) 
dx 


K 1 

$n ~ 2 

K 1 


df dK ' 

— - l-K-’f 

dx dx 


9{t)+ 


dR 2 _ dC 
— nC 

dx dx 


*n’C^+ 

dx 


(4.19) 


dK 


$n 


-2 


dM 


q- *n 


_ 3 ^q 

dx 


dx dx 

The key to the effectiveness of this mode acceleration sensitivity method is the usage 
of the exact K~ 1 in the calculation of the K 1 dK/dx$f) 2 and K 1 dM/dx# 
terms. The explicit calculation of these terms expands the basis beyond the span 
of the modes in a manner similar to the pseudo-static term in the modified modal 
method described in the previous section. 


Using equation (4.19) the stress sensitivities can be calculated as 

dcr „ du dS 

— = S 1 u 

dx dx dx 

where u is obtained from equation (2.11). 

It is worthwhile to contrast the sensitivity approach of equations (4.19) with 
an alternate approach of a fixed mode, overall forward difference method with the 
response quantities calculated using a mode acceleration method. This overall for- 
ward difference approach has one obvious drawback. The mode acceleration method 
requires the costly factorization of K for the perturbed design. So, much of the cost 
savings achieved by keeping the modes fixed is lost. A second defect of the overall 
forward difference method is that it is not as effective as the method of equations 
(4.19). The mode acceleration method in the overall finite difference procedure 
provides a good approximation to the first term (pseudo-static term) in equations 
(4.19). However, the key effects of the K _1 dK/dx3>0~ 2 and the K _1 dM/dx# 
terms are completely neglected. 


(4.20) 
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Chapter 5 


Numerical Studies 

The different transient response methods described in Chapter 2 and the sen- 
sitivity calculation methods described in Chapter 4 are applied to three example 
problems in this chapter. The three example problems are small but they all have 
certain characteristics which complicate the dynamics and sensitivity calculations. 
The first example, the five- span beam has relatively closely spaced frequencies and 
is loaded with a moment applied at a single point. As a result, many modes partic- 
ipate in the dynamic response. The second example, the delta wing, is loaded with 
a uniform pressure load. Although the higher frequency modes are not significantly 
excited by this loading, the analysis is complicated by the laminated plate elements 
in the model and the sensitivity analysis is complicated by the lamina thickness 
design variables considered. The third example is a cantilever beam with a stepped 
cross section loaded with an applied rotation at the root. This loading is inertial, 
depends on the mass, and therefore also depends on the values of the design vari- 
ables. The first two examples consider point mass and standard thickness design 
variables. The cantilever beam example also includes so-called “shape” design vari- 
ables (section lengths) that are known to cause difficulties in the sensitivity analysis 
in some cases. 

One of the key questions addressed in this chapter is how well a particular set 
of basis vectors represents the full system of order n g . This full system, however, 
is the result of a particular finite element discretization. Thus the accuracy of 
the response or sensitivities as a function of the finite element mesh is also an 
appropriate question. This question is especially important when a large number 
of basis vectors (n r close to n g ) are required for an accurate solution in a problem 
with a given finite element mesh. Either the basis vectors are doing a very poor 
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job of spanning the solution space or the loading is legitimately exciting this high 
frequency behavior. 

In this chapter two terms will be used to describe these studies which consider 
the dynamic response as a function of the number of basis vectors or the number 
of finite elements in the model. The effect of the number of basis vectors on the 
accuracy of the response or sensitivities for a given finite element mesh is called 
a “modal convergence” study. In this case, the goal is for the n r basis vectors to 
provide an accurate solution to the approximate equations of order n g . The question 
of whether the finite element model associated with this system is an accurate 
representation of the continuum is addressed in a “mesh convergence” study. In 
some cases it will be shown that the modal convergence is strongly related to the 
mesh convergence. That is, when a large number of basis vectors are required for an 
accurate solution for a given finite element mesh, the finite element mesh is doing 
a poor job of representing the continuum. In other cases, even though the finite 
element mesh is providing an accurate representation of the continuum, some sets 
of basis vectors are doing a poor job of representing the response or sensitivities for 
this n 9 th order system. 

Several additional comments on the concept of a “modal convergence” study are 
in order. Clearly the use of the term convergence is imprecise because the accuracy 
of the approximate solution with different numbers of modes is compared only 
against the finite degree-of-freedom solution rather than the continuum solution. 
However, it is assumed that an “acceptable” finite element model must do a good job 
of representing the low frequency modes of the structure. Therefore, the accuracy 
of the dynamic solution with a small number of modes from the finite degree-of- 
freedom model is a very reasonable approximation to the accuracy of the dynamic 
solution with a small number of modes calculated from a continuum model. Thus 
the convergence of the solution as a function of the number of modes calculated 
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from the finite element model is a reasonable approximation to the true modal 
convergence obtained when the modes are calculated from the continuum model as 
long as the number of modes considered is small. Furthermore, if the number of 
modes required for an accurate calculation of either the response or sensitivities is 
not small, the basis vectors or the method will be considered poor. 

5.1. Five-span beam Example 

The first example considered is a five-span, planar beam example taken from 
reference [33] and shown in figure 5.1. An initial investigation of the displacement 
transient response of this problem was also considered in [34]. In most of the stud- 
ies, the beam is modeled with three beam finite elements per span resulting in 26 
unconstrained degrees-of-freedom. The efFect of finite element discretization is con- 
sidered by developing alternate models with 6, 9, 12, etc. elements per span. As 
shown in figure 5.1, translational and rotational viscous dampers were also added 
to the beam. These devices are representative of velocity feedback controllers which 
might be added to flexible structures. Cases with and without dampers were con- 
sidered. The numerical values of the damping coefficients from ref. [33] of c\ — .008 
sec-lbf/in and ci — 1.2 sec-lbf were used. In one example, modal damping with 
= .005, which is intended represent typical structural damping, was used instead 
of the discrete dampers. A case was also considered where a 1.0 lb mass (approx- 
imately 20% of the beam’s mass) was added to the beam at the location of the 
translational damper. The eigenvalues for three cases using the three element per 
span model are shown in table 5.1. The additional point mass has a significant 
efFect on the frequencies while the dampers have little effect. The effect on frequen- 
cies of increasing the number of elements per span in the finite element model is 
shown in table 5.2. It can be seen that the lowest ten frequencies are fairly well 
converged even for the model with 3 elements per span. In the transient analysis, 
the applied loading for all problems consisted of a point moment of .04405 in-lbs 
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applied at the right end of the beam. Two different time functions for this load, a 
step and a ramp, (shown in figure 5.1) were considered. 


u u u u u 

A AcJaJ/ a-4/ ) m 


Beam cross section 
H-2.0 ln.~H j 


|*30 in=*| 


M .004405 
In.-lbs 


Ramp load 


, M .004405 
in.-lbs 


.2 


t, secs 


E = 10 x 10 5 lb/in? 


Step load 


h = .0625 In. 


t, secs 


p = .28 Ib/in. 

c, = .008 sec - Ib/in. 

= 1.2 sec - lb 

Figure 5.1. Five-span beam with applied end moment. 


Table 5.1. Eigenvalues For Three Five-span beam Cases 


Mode 

Undamped 

Damped With Point Dampers 

Undamped 
With Point Mass 


Frequency, Hz. 

Frequency, Hz. 

Damping Ratio 

Frequency, Hz. 

■B 

1.1707 

1.2210 

.0851 

.9401 


1.2991 

1.2926 

.0352 

1.2594 


1.6254 

1.6298 


1.5445 


2.0491 

2.0910 

.0590 

1.8005 

5 

2.4628 

2.5497 

.0958 

2.3729 

6 

4.7343 

4.8426 

.0044 

4.2327 

mm 

5.0105 

4.9785 

.0413 

4.8858 


5.6472 

5.7703 


5.6400 

H 

6.4153 

6.4178 


5.9261 

n 

7.1274 

7.2229 

.0193 

6.8762 
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Table 5.2. Beam Frequencies With Different Numbers of Finite Elements Per Span 


Mode 

Frequency, Hz 

3-Elements 

6-Elements 

9-Elements 

12-Elements 

1 

1.1707 

1.1698 

1.1698 

1.1698 

2 

1.2991 

1.2979 

1.2978 

1.2978 

3 

1.6254 

1.6230 

1.6229 

1.6229 

4 

2.0491 

2.0445 

2.0442 

2.0442 

5 

2.4628 

2.4547 

2.4542 

2.4542 

6 

4.7343 

4.6828 

4.6798 

4.6792 

7 

5.0105 

4.9504 

4.9469 

4.9462 

8 

5.6472 

5.5652 

5.5601 

5.5593 

9 

6.4153 

6.3053 

6.2980 

6.2967 

10 

7.1274 

6.9974 

6.9874 

6.9857 


5.1.1. Beam Dynamic Response 

The first part of this study focused on the transient response of the beam using 
the mode displacement, mode acceleration, static mode, and Ritz-Wilson-Lanczos 
(RWL) methods. Displacement, velocity, acceleration, and stress resultant response 
quantities are considered. For this beam example all of these response quantities 
are taken at a location 10.0 inches from the left end of each span. This point is the 
end of the first element in each span when three elements per span are used in the 
model. 

5. 1.1.1. Character of the response 

In the first case the ramp loading was applied to the undamped beam modeled 
with three elements per span. All 26 modes were used in the analysis. Time histories 
of selected displacement, velocity, acceleration, and bending moment components 
are shown in figures 5.2, 5.3, 5.4, and 5.5, respectively. The displacement history 
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(fig. 5.2) is relatively smooth indicating that only the low frequency modes of the 
beam are contributing to the response. The velocity and bending moment response 
histories are more jagged indicating participation by higher frequency modes. The 
acceleration history (fig. 5.4) is extremely jagged with contributions from the high- 
est frequency modes represented by the finite element model. 



Time, secs 

Figure 5.2. Time history of displacement u 2 for five-span beam subjected to tran- 
sient end moment (Ramp load, undamped beam). 


The impulsive nature of the step load makes the higher frequency modes much 
more important. This can be seen in figure 5.6 where the time history of velocity 
in the second span (see fig. 5.1), U 2 is shown. By comparing this velocity history 
with that in fig. 5.3 the increased importance of the high frequency modes becomes 
obvious. 

The addition of the point dampers shown in fig. 5.1, on the other hand, tends 
to reduce the importance of the high frequency modes. This is shown in fig. 5.7 
where again U 2 is shown. Comparing this to fig. 5.3 it can be seen that the velocity 
history for the damped case is significantly smoother than for the undamped case. 
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Time, secs 

Figure 5.3. Time history of velocity U2 for five-span beam subjected to transient 
end moment (Ramp load, undamped beam). 



.00 .50 1 .00 1 .50 2.00 2.50 3.00 


Time, secs 

Figure 5.4. Time history of acceleration u 2 for five-span beam subjected to tran- 
sient end moment (Ramp load, undamped beam). 
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Time, secs 

Figure 5.5. Time history of bending moment in span five for five-span beam sub- 
jected to transient end moment (Ramp load, undamped beam). 


u 2 

in. /sec 



Figure 5.6. 


Time, secs 

Time history of velocity U 2 for five-span beam subjected to transient 
end moment (Step load, undamped beam). 
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This changing character of the time histories with temporal or spatial differ- 
entiation of the response function or the addition of dampers is expected. The 
implications of this phenomenon on calculating the sensitivities of these response 
quantities will be discussed below. 



Time, secs 

Figure 5.7. Time history of velocity u? for five-span beam subjected to transient 
end moment (Ramp load, damped beam). 
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5. 1.1. 2. Modal convergence 


When vibration inodes or other functions are used to reduce the basis in a 
transient response problem (eq. (2.4)), the key question is how many modes are 
required for an accurate solution. This section addresses that question for the 
five-span beam example with the response calculated using mode displacement, 
mode acceleration, static mode, and RWL methods. Unless otherwise stated, all of 
the response quantities are considered at critical times selected using the methods 
discussed in Chapter 3. 

The baseline case of ramp loading applied to the undamped beam modeled 
with three elements per span is considered first. Figure 5.8 shows the convergence 
of selected displacements at critical points as a function of number of modes. The 
displacement/critical point combinations were selected to be representative of both 
the largest and smallest critical values. In figure 5.8 and in all the other figures 
showing convergence of response quantities or sensitivities, the figure key indicates 
the quantity followed by the time of occurrence in seconds. In all cases the con- 
vergence is very good with approximately ten modes yielding a converged solution. 
Figure 5.9 show a similar plot for velocities. Again the convergence is good. The 
modal convergence for accelerations, however, is poor as shown in fig. 5.10. Figures 
5.11 and 5.12 show the modal convergence of selected bending moments and shear 
forces respectively. Again the convergence is poor. 

To possibly alleviate this poor convergence, the alternate reduction methods 
discussed in Chapter 2 were applied to this problem. The modal convergence for 
displacements calculated using the mode acceleration method (fig. 5.13) is even 
better than that found using the mode displacement method. The convergence 
of bending moments and shear forces has improved dramatically from the mode 
displacement results as can be seen in figs. 5.14 and 5.15. As mentioned in chapter 
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Figure 5.8. Modal convergence of selected displacements for the five-span beam 
(Ramp load, undamped beam, mode displacement method). 
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Figure 5.9. Modal convergence of selected velocities for the five-span beam (Ramp 
load, undamped beam, mode displacement method). 
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Figure 5.10. Modal convergence of selected accelerations for the five-span beam 
(Ramp load, undamped beam, mode displacement method). 
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Figure 5.11. Modal convergence of selected bending moments for the five-span 
beam (Ramp load, undamped beam, mode displacement method). 
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A similar improvement was noted using the static mode method. As an example 
consider the excellent convergence of shear forces shown in fig. 5.16. However, 
the addition of the static solution provides no improvement in the convergence of 
acceleration as shown in fig. 5.17. 

The RWL method is attractive because of the significantly reduced cost of 
calculating the vectors compared with solving the eigenproblem. In this five-span 
beam example the modal convergence is also as good as the mode acceleration or 
static mode methods. The good convergence of the shear forces is shown in fig. 
5.18. Like the other reduction methods, however, the convergence of accelerations 
is poor (fig. 5.19). 
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The modal convergence of the response quantities for the step loaded case is 
generally much poorer than for the ramp loaded case. The convergence of the 
displacements is reasonably good. Convergence of velocities, accelerations, and 
stresses, however, is poor. This poor convergence is not surprising considering the 
“jaggedness” of the velocity time history shown in figure 5.6. As an example, two 
figures plotting the convergence of bending moments as a function of the number 
of modes are shown. The first, figure 5.20, shows the bending moments calculated 
using the mode displacement method. Convergence is poor but this is not surprising 
since the convergence was poor with the mode displacement method for the ramp 
loaded case (figure 5.11). For this ramp loaded case, the convergence of the bending 
moments improved dramatically when the mode acceleration method was used as 
can be seen in figure 5.14. Although convergence is improved for the step loaded 
case by using the mode acceleration method (figure 5.21), the convergence is still 
fairly poor. 
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Figure 5.13. Modal convergence of selected displacements for the five-span beam 
(Ramp load, undamped beam, mode acceleration method). 
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Figure 5.14. Modal convergence of selected bending moments for the five-span 
beam (Ramp load, undamped beam, mode acceleration method). 
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Figure 5.15. Modal convergence of selected shear forces for the five-span beam 
(Ramp load, undamped beam, mode acceleration method). 
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Figure 5.16. Modal convergence of selected shear forces for the five-span beam 
(Ramp load, undamped beam, static mode method). 
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Figure 5.17. Modal convergence of selected accelerations for the five-span beam 
(Ramp load, undamped beam, static mode method). 
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Figure 5.18. Modal convergence of selected shear forces for the five-span beam 
(Ramp load, undamped beam, Ritz-Wilson-Lanczos method). 



(1.25) 

(1.58) 

(2.30) 

(2.84) 


Number of modes 

Figure 5.19. Modal convergence of selected accelerations for the five-span beam 
(Ramp load, undamped beam, Ritz-Wilson-Lanczos method). 
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Figure 5.20. Modal convergence of selected bending moments for the five-span 
beam (Step load, undamped beam, mode displacement method). 
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Figure 5.21. Modal convergence of selected bending moments for the five-span 
beam (Step load, undamped beam, mode acceleration method). 
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Judging from the velocity time history in figure 5.7, it might be expected that 
including damping would improve the modal convergence of the response quanti- 
ties. For the ramp loaded, undamped case, the poorest convergence was for the 
accelerations (figure 5.10). For the ramp loaded case with discrete dampers there is 
an improvement in modal convergence as seen in figure 5.22. For the case with .5% 
modal damping there is also a slight improvement in modal convergence. However, 
in neither case does the damping completely alleviate the poor convergence. 



Number of modes 

Figure 5.22. Modal convergence of selected accelerations for the five-span beam 
(Ramp load, discretely damped beam, mode displacement method). 


All of the above convergence results are at critical points located by the method 
described in Chapter 3. When a different number of modes is used in the analysis, 
the critical time for a particular critical point usually shifts slightly. Consequently, 
the results for a given response quantity/critical point combination occur at different 
times depending on the number of modes used in the analysis (the values shown 
in parenthesis in the figures are for the most refined solution). It is natural to ask 
whether response quantities at fixed times plotted as a function of number of modes 
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would show similar convergence. Figures 5.23 and 5.24 show the modal convergence 
of selected velocities and bending moments, respectively, at fixed times. The mode 
displacement method was used in the analyses. The particular response quantities 
and times were selected to span the range between largest positive and negative 
values. As can be seen in fig. 5.23, the convergence of velocities is good. From 
fig. 5.24 it can be seen that the convergence of the bending moments is poor and 
remarkably similar to the critical point convergence results (fig. 5.11). Thus it 
would appear that the critical point constraint formulation does not significantly 
affect the modal convergence of the response. 
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Figure 5.23. Modal convergence of selected velocities for the five-span beam at 
fixed time points (Ramp load, undamped beam, mode displacement 
method). 
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Figure 5.24. Modal convergence of selected bending moments for the five-span 
beam at fixed time points (Ramp load, undamped beam, mode dis- 
placement method). 


5. 1.1. 3. Mesh convergence 


Table 5.2 shows the convergence of the lowest ten frequencies as a function of 
the number of elements used to model each span of the beam. The convergence of 
these lower frequencies is rapid. The convergence of various response quantities as 
a function of the mesh is also a concern. 

The modal convergence can not be uncoupled from the mesh convergence. 
This was discussed in ref. [33] for the derivatives of damping ratios calculated 
using undamped vibration modes. For several cases, the modal convergence of the 
derivatives was poor. As the mesh was refined, convergence was achieved only when 
almost all of the available modes were used in calculating the damping ratio. Clearly, 
this is an example where the modal basis provides a very poor approximation to 
the actual solution. 
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Figure 5.25 shows the modal convergence of shear force for the five-span beam 
modeled with six elements per span and the transient analysis performed using the 
mode displacement method. The convergence for this case is just as poor as for the 
three-element-per-span case shown in fig. 5.12. However, a plot of convergence of 
this shear forces as a function of number of elements per span, when all modes are 
used in each analysis, (fig. 5.26) shows good convergence. Clearly, the convergence 
of shear forces for the ramp-loaded five-span beam is similar to that reported for 
derivatives of damping ratios in ref. [33]; the vibration modes are simply doing a 
poor job of representing the solution. 
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Figure 5.25. Modal convergence of selected shear forces for the five-span beam 
modeled with six elements per span (Ramp load, undamped beam, 
mode displacement method). 


Figure 5.27 which shows the convergence of accelerations for the step loaded 
beam as a function of elements per span indicates a different behavior, however. 
Here, the very poor convergence is due to the higher frequency modes being excited 
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Figure 5.26. Convergence of selected shear forces as a function of number of finite 
elements per span for the five-span beam (Ramp load, undamped 
beam, mode displacement method). 


by the step loading. As the mesh is refined the number of high frequency modes 
increases and these continue to have a significant contribution to the acceleration. 


In evaluating the accuracy of the sensitivity calculation procedures in the next 
section, particular attention must be paid to the convergence characteristics. Some 
convergence problems such as those caused by the use of vibration mode shapes can 
be improved by the use of alternate basis functions. However, other convergence 
problems, such as for the accelerations in the step-loaded case are inherent in the 
problem definition. 
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Figure 5.27. Convergence of selected accelerations as a function of number of finite 
elements per span for the five-span beam (Step load, undamped beam, 
mode displacement method). 


5.1.2. Sensitivities of Beam Dynamic Response 


In the previous section the transient response of the five-span beam was con- 
sidered in detail. In this section the calculation of sensitivities of displacements, 
velocities, accelerations, and stresses with respect to various design variables is 
considered. 

5. 1.2.1. Design variables 


Two different classes of design variables were considered. The first design vari- 
able is a concentrated mass, m (initially zero) at the location of the translational 
damper. This design variable was also considered in [33]. The derivatives of the 
system mass and stiffness matrices with respect to this design variable are constant 
and zero, respectively. As a consequence, the derivatives of the system matrices re- 
quired in the semi-analytical methods can be calculated exactly by a simple forward 
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difference operator. The beam thicknesses in each of the five spans were also design 
variables. Derivatives with respect to the five thickness design variables showed 
similar characteristics. Herein, results for derivatives with respect to the thickness 
in the rightmost span, / 15 , along with derivatives with respect to the point mass, 
m, are presented. 

5. 1 . 2 . 2. Effect of finite difference step size 

The methods described in Chapter 4 for calculating sensitivities all rely on 
finite difference operators at some stage in the algorithm. The forward and central 
difference methods rely on the operators in equations (4.1) and (4.2) to calculate 
derivatives of response quantities. In the semi-analytical methods the derivatives of 
the system matrices are calculated using the forward difference operator in equation 
(4.1). In all cases the finite difference step size must be selected so that the operator 
provides a reasonable approximation to the derivative. If the step size is too large, 
the error due to truncating the series approximation of the derivative is large. If 
the step size is too small, the numerical condition error in performing the function 
evaluations (dynamic analyses) becomes large. 

To assess the effect of step size on the calculation of sensitivities for the five- 
span beam, derivatives were calculated using the three methods with various step 
sizes. In this study the beam was undamped and the ramp loading was applied. 
All 26 vibration modes were included in the analysis. Figure 5.28 shows the esti- 
mated derivative of displacement u 2 at critical point number 5 with respect to the 
point mass design variable, m, as a function of step size. As mentioned above, the 
derivatives of the system matrices with respect to this design variable can be calcu- 
lated exactly in the semi-analytical method. As a result, the derivative estimated 
using the semi- analytical method is approximately constant for the six-order-of- 
magnitude change in step size shown in the figure. The central difference method 
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uses the higher order operator and provides good accuracy over most of the step size 
range shown in the figure. The forward difference operator provides good accuracy 
with the smaller step sizes but begins to diverge earlier for the larger step sizes than 
the central difference method. 


du2 

dm 



O Central difference 
□ Semi-analytical 
A Forward difference 


Finite difference step size, lbs 

Figure 5.28. Effect of finite difference step size on the accuracy of a displacement 
derivative approximation with respect to the point mass design vari- 
able (Ramp load, undamped beam). 


Figure 5.29 shows the estimated derivative of displacement u\ at critical point 
number 5 with respect to the right-most span thickness, h 5 , as a function of step size. 
In this case the system mass matrix is a linear function of this design variable and 
its derivative can be represented exactly by the forward difference operator. The 
system stiffness matrix is a cubic function of this design variable and its derivative 
can only be approximated by the forward difference operator. Still, the derivative 
approximation computed by the semi-analytical method is very accurate except for 
the largest step size and is no worse for this case than the much more costly central 
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O Central difference 
□ Semi-analytical 
A Forward difference 


Finite difference step size, in. 

Figure 5.29. Effect of finite difference step size on the accuracy of a displacement 
derivative approximation with respect to the thickness design variable 
(Ramp load, undamped beam). 


difference method. Again, the forward difference operator results in substantial 
errors for the larger step sizes. 


Because this example has a relatively small number of degrees of freedom there 
is little condition error when small step sizes are used. To assess the effects of 
condition error which would occur for larger problems, the derivative approxima- 
tions for the five-span beam problem were also calculated using 32- bit floating point 
precision compared with the 60-bit precision used in the studies described above. 
The estimated derivative of displacement U 2 at critical point number 5 with respect 
to the point mass is plotted as a function of finite difference step size in figure 
5.30. Derivative approximations are calculated using the semi-analytical method, 
the central difference method, and the forward difference method. For the larger 
step sizes, the results from all three methods agree well with those calculated using 
60-bit precision. For step sizes smaller than 10 -4 there is considerable error in the 
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Figure 5.30. Effect of finite difference step size on the accuracy of a displacement 
derivative with respect to the point mass design variable. Calcu- 
lations performed using 32 bit precision. (Ramp load, undamped 
beam). 

derivative approximations calculated using the forward and central difference meth- 
ods. The derivative approximations calculated using the semi- analytical method, 
however, are highly accurate over the entire range of step sizes shown in the figure. 


5. 1.2. 3. Modcil convergence of sensitivities 


The first case considered is the undamped beam with the ramp load. Figure 
5.31 shows the convergence of selected estimated derivatives of displacements with 
respect to m at various critical points. The mode displacement method was used and 
the derivative approximations were calculated using the central difference operator 
with updated modes. The convergence is good although slightly poorer than the 
convergence of the displacements themselves (fig. 5.8). The convergence of the 
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Number of modes 

Figure 5.31. Modal convergence of derivatives of selected displacements with re- 
spect to the mass design variable (Ramp load, undamped beam, mode 
displacement method, central difference operator). 

estimated displacement derivatives with respect to the thickness design variable is 
similar. 

Although the modal convergence of the velocities for this case is good (fig. 5.9), 
the convergence of selected estimated derivatives of velocity with respect to m is 
generally poor (fig. 5.32). Considering the poor convergence of the accelerations 
shown in fig. 5.10 it is not surprising that the convergence of the sensitivities of 
the accelerations is also very poor. From fig. 5.33 it can be seen that the derivative 
approximations of the four selected critical point accelerations with respect to the 
thickness design variable are essentially not converging with increasing number of 
modes. It should be pointed out again that these derivative approximations of 
velocity and acceleration are calculated using the central difference method and 
updated mode shapes thus minimizing the numerical errors. The poor convergence 
exhibited in figs. 5.32 and 5.33 is due to the poor approximation of the sensitivities 
by the mode shapes. 
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Figure 5.32. Modal convergence of derivative approximations of selected veloci- 
ties with respect to the mass design variable (Ramp load, undamped 
beam, mode displacement method, central difference operator). 

Similar modal convergence behavior is observed for sensitivities of the stress 
resultants. This is consistent with the poor convergence of the stress resultants 
calculated using the mode displacement method (figs. 5.11 and 5.12). Figure 5.34 
shows the poor convergence of derivative approximations of selected bending mo- 
ments with respect to the thickness design variable. It can be seen that the conver- 
gence of the bending moment derivative approximation in the rightmost span with 
respect to the thickness in the rightmost span (dM& / dh{) is especially poor. 

It was shown in the previous section that several approaches are available for 
overcoming the poor convergence of bending moments and shear forces in this beam 
example. The mode acceleration, static mode, and RWL methods all produced good 
modal convergence of bending moments and shears as shown in figures 5.14, 5.15, 
5.16, and 5.18. Unfortunately this type of dramatic improvement does not occur 
for the sensitivities of the stress resultants. Figure 5.35 shows the convergence 
of the bending moment derivative approximations with respect to the thickness 
design variable where the analysis was done using the RWL method. As in the 
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Figure 5.33. Modal convergence of derivative approximations of selected accelera- 
tions with respect to the thickness design variable (Ramp load, un- 
damped beam, mode displacement method, central difference opera- 
tor). 


studies discussed above using the mode displacement method, the sensitivities were 
calculated using the central difference operator with the basis vectors updated for 
the perturbed design. Convergence of dM^/dhs is somewhat improved compared 
to the mode displacement case. Other quantities show convergence similar to the 
mode displacement case; none of these convergence histories can be described as 
good. Convergence of the shear force derivative approximations using the RWL 
method is considerably worse than for the bending moments. 


The semi-analytical methods have also been used for calculating sensitivities of 
stress resultants. Figure 5.36 shows the convergence of bending moment derivative 
approximations with respect to the mass design variable calculated using the fixed- 
mode, semi-analytical method and RWL vectors. The convergence is very similar, 
especially for larger numbers of modes, to that of the central difference method. 
The mode acceleration, semi-analytical method and the semi- analytical method 
with approximate d&/dx were also tried. Again, the modal convergence curves had 
the same “jaggedness” as for previous cases. 
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Figure 5.34. Modal convergence of derivative approximations of selected bend- 
ing moments with respect to the thickness design variable (Ramp 
load, undamped beam, mode displacement method, central difference 
operator). 



Number of modes 

Figure 5.35. Modal convergence of derivative approximations of selected bending 
moments with respect to the thickness design variable (Ramp load, 
undamped beam, RWL method, central difference operator). 
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Figure 5.36. Modal convergence of derivative approximations of selected bending 
moments with respect to the mass design variable (Ramp load, un- 
damped beam, RWL method, semi-analytical formulation). 


Considering the above difficulties with modal convergence for the ramp loaded 
cases, especially poor convergence would be expected for the step loaded case. For 
the ramp loaded case the convergence of displacement derivative approximations 
with respect to the mass design variable was reasonably good (figure 5.31). For 
the step loaded case the modal convergence of the same displacement sensitivities 
is poor as shown in figure 5.37. The modal convergence of higher order sensitivities 
(velocities, accelerations, and stresses) is extremely poor. 


Adding damping slightly improved the modal convergence of the response quan- 
tities but did not completely alleviate the convergence problem. The result for 
sensitivities is similar. Figure 5.38 shows the convergence of velocity sensitivities 
for the discretely damped case. The convergence is slightly improved over the un- 
damped case shown in figure 5.32 but the curves are still fairly jagged. Convergence 
of sensitivities for the case with modal damping is also very similar to that in the 
undamped case. 
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Figure 5.37. Modal convergence of derivative approximations of selected displace- 
ments with respect to the mass design variable (Step load, undamped 
beam, mode acceleration method, semi-analytical method). 
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Figure 5.38. Modal convergence of derivative approximations of selected veloci- 
ties with respect to the mass design variable (Ramp load, discretely 
damped beam, mode acceleration method, semi-analytical method). 
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5. 1.2. 4. Mesh convergence of sensitivities 


Just as for the response quantities, additional insight can be obtained by con- 
sidering the convergence of their sensitivities with increasing number of elements 
per bay. The case of the stress resultants will be considered since they were shown 
to converge well with mesh refinement (fig. 5.26) but poorly as a function of num- 
ber of vibration modes used in the analysis. Figure 5.39 shows the convergence of 
shear force derivative approximations with respect to the mass design variable as a 
function of number of elements per bay. Surprisingly, the convergence is extremely 
poor. 
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Figure 5.39. Convergence of derivative approximations of selected shear forces 
with respect to the mass design variable as a function of number of 
elements per span (Ramp load, undamped beam, mode displacement 
method, central difference operator). 
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5. 1 . 2 . 5. Fixed versus updated modes in sensitivity calculations 


As mentioned above, the computational cost of updating the vibration modes 
for the perturbed analyses is substantial. The question of whether the modes from 
an initial design can be used in a finite-difference-based procedure to calculate sensi- 
tivities of the transient behavior has received considerable attention in the literature. 
In ref. [35] it was shown that there is a substantial difference in the derivatives of 
aircraft flutter speeds when fixed modes are used rather than the updated modes. 
In ref. [33], however, there was little difference in the derivatives of damping ratios 
for the five-span beam when either fixed or updated modes were used. This was 
investigated here using the same five-span, undamped beam under the step load. As 
shown in figure 5.37 where the derivative approximations were calculated using the 
semi-analytical, mode acceleration method, convergence with respect to the number 
of modes is very slow. Figure 5.40 shows the modal convergence of derivative ap- 
proximations of selected displacements with respect to /15 calculated using forward 
difference procedures. Results with both fixed and updated vibration modes are 
shown. Again, the convergence as a function of number of modes is poor. However, 
for till three derivative approximations, the results are nearly the same for both the 
fixed mode and updated mode cases. 
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Figure 5.40. Modal convergence of derivative approximations of selected displace- 
ments with respect to the mass design variable calculated with both 
fixed and updated vibration modes (Step load, undamped beam, for- 
ward difference methods). 
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5.2. Composite Delta wing Example 


The second example considered is an aircraft delta wing with laminated com- 
posite cover skins taken from ref. [36] and described in detail in ref. [37]. The finite 
element model of this structure is shown in fig. 5.41. Since the wing is geometrically 
symmetric about the midplane through its thickness (X-Y plane), only the upper 
half of the wing is modeled and boundary conditions enforcing anti-symmetric mo- 
tion are imposed on the joints lying in the X-Y plane. The wing is also cantilevered 
at the root. The model contains a total of 88 joints with a total of 140 unconstrained 
degrees-of-freedom. The webs in the wing are made of titanium and are modeled 
with 70 shear panel finite elements along with rod elements through the thickness 
of the wing. The cover skins are made of a moderate-modulus ( E \ = 21 x 10 6 psi) 
graphite-epoxy material with 0°, ±45°, and 90° lamina where the 0° material runs 
spanwise along the wing. These cover skins are modeled with membrane finite ele- 
ments so only the total thicknesses (and not the stacking sequence of plies) of each 
la mi na are important. The structural mass is 6003 lbs but most of the wing mass is 
due to a fuel mass of 93650 lbs distributed over the joints. The spatial distribution 
of the load is the same as the static load from ref. [37] and is roughly equivalent to 
a 144 psf pressure load on the wing skin. A step loading function was used as the 
time function for all cases. The lowest ten vibration frequencies for the wing are 
shown in table 5.3. Damping is accounted for by assuming .5% of critical damping 
for all modes. 

5.2.1. Wing Dynamic Response 

The character of the dynamic response of the delta wing is considerably dif- 
ferent than that of the five-span beam. Shown in figure 5.42 is a time history of 
acceleration at the wing tip. Although 64 modes were included in the analysis, it is 
evident from figure 5.42 that only the low frequency modes are being excited. The 
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Figure 5.41. Finite element model of composite delta wing. 

Table 5.3. Lowest 10 Vibration Frequencies For the Delta Wing 


— 

Mode 

Frequency, Hz 

1 

2.055 

2 

2.765 

3 

4.104 

4 

4.913 

5 

5.920 

6 

6.944 

7 

7.451 

8 

8.421 

9 

9.583 

10 

9.880 


same is true for stresses as shown in the time history of figure 5.43. Shown in fig. 
5.43 is which is a typical shear stress in a web. As can be seen, there is a small 
amount of higher frequency response superimposed on the predominant response 
frequency. However, the time history exhibits none of the high-frequency response 
present in the five-span beam. 
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5. 2. 1.1. Modal convergence 


In contrast to the five-span beam example, the modal convergence of all re- 
sponse quantities considered for the delta wing is quite good. Shown in figures 
5.44 and 5.45 are modal convergence plots for selected accelerations and stresses 
at critical points calculated using the mode displacement method. A converged 
solution is reached with approximately twenty or less modes for all of the response 
quantities shown. Convergence is also good for response quantities when the mode 
acceleration or RWL methods are used instead of the mode displacement method. 
Shown in figure 5.46 is a convergence plot for the same stresses shown in figure 5.45 
but calculated using RWL vectors. 
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Figure 5.44. Modal convergence of tip accelerations for the delta wing. 
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Figure 5.45. Modal convergence of selected stresses for the delta wing calculated 
using the mode displacement method. 
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Figure 5.46. Modal convergence of selected stresses for the delta wing calculated 
using the RWL method. 
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5.2.2. Sensitivities of Wing Dynamic Response 
5 . 2 . 2 . 1 . Design variables 

The design variable definitions are the same as those in ref. [37] and are shown 
in figure 5.47. As can be seen in figure 5.47 the skin is broken up into 16 regions. 
In each region there are three design variables — the thickness of the 0° lamina, the 
thickness of the 90° lamina, and the thickness of the ±45° lamina. These design 
variables will be denoted t l 0 where i denotes the region of the wing skin, and 9 is 
either 0, 90, or 45 depending on the lamina orientation. Also shown in figure 5.47 
are the 12 design variables controlling the thickness of the webs. These will be 
denoted t l w where i denotes the particular web region. 



Figure 5.47. Design variable definitions for the delta wing example. 
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In calculating sensitivities of various response quantities, only a small subset 
of these design variables were considered. Specifically, derivatives of selected dis- 
placement, velocity, acceleration, and stress quantities were calculated with respect 

t 0 y2 ,13 ,16 ,16 ,6 _ n J ,10 
lO Iq, Iqq, Iq J l 90» L wt aIlu 1 W ' 

5. 2. 2. 2. Effect of finite difference step size 

Compared to the five-span beam example, the system matrices for the delta 
wing are larger and have a more complicated connectivity. Since many of the signifi- 
cant operations in the transient response analysis operate directly on these matrices, 
there is considerable potential for accumulating roundoff error. This roundoff error 
along with the truncation error in the finite difference expressions is a concern in 
selecting a step size for a finite difference approximation to a derivative. 

A study was performed to consider the effect of step size in the forward differ- 
ence and central difference methods for the delta wing. Figure 5.48 shows derivative 
approximations of the wing tip acceleration at critical points with respect to selected 
thickness design variables as a function of the finite difference step size used. As 
seen in the figure the step size was varied by factors of ten from 10 -7 to 10 -2 . The 
central difference method was not used with the 10 -2 step size because the backward 
perturbation from the nominal design would result in negative member thickness. 
One significant observation is that the acceptable step size range for the forward 
difference method is small — approximately two decades. A second observation is 
that the behavior of the central difference method as a function of step size is sur- 
prisingly good. It is expected that, for larger step sizes (10 -3 ), the central difference 
method results in less error than the forward difference method. The unexpected, 
superior performance of the central difference method for the smaller step sizes is 
probably due to the symmetry of the difference operator. The roundoff errors that 
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Figure 5.48. Effect of finite difference step size on the accuracy of tip displacement 
derivative approximations calculated using the forward and central 
difference methods with fixed modes (Delta wing example). 


occur with the positive and negative perturbations tend to cancel each other and 
thus produce the better-than-expected, accurate values for the sensitivities. 


The situation is similar for selected stress sensitivities shown in figure 5.49. 
Most of the curves for the forward difference case have a small acceptable step 
size range. This is especially obvious for the derivative dcr\\/dt \ ® where 10 -5 is 
the apparent choice for step size. It should also be mentioned here that these 
calculations were performed using 64 bit arithmetic. In the five-span beam example, 
the effect of step size on displacement derivatives was not as severe even though, 
for one case, these were calculated with predominantly 32 bit arithmetic (fig. 5.30). 

The simple approach of selecting a single step size for use with all response 
quantities and all design variables was used here. This approach has the obvious 
advantage of simplicity but very questionable validity for the forward difference 
method and this delta wing example. From figure 5.49 there is significant error in 
d^li/dtlg if greater than 10 -4 is used as the step size. However, if less than 10 -5 
is used instead, dutip/dvft o is in error. 
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Figure 5.49. 

As noted above, the central difference method improves the range of accept- 
able step sizes but at the cost of an additional analysis for each design variable. 
Alternatively, the semi-analytical method is particularly attractive for this delta 
wing example. The stiffness matrices of the membrane and shear panel finite ele- 
ments are linear functions of the thickness design variables. Thus large values of 
the step size can be used to effectively eliminate the roundoff errors in generating 
the derivative approximations of the stiffness and mass matrices required for the 
semi-analytical method. 

5. 2. 2. 3. Modal convergence of sensitivities 

Unlike the five-span beam example, the modal convergence of the displacement, 
velocity, and acceleration derivatives for the delta wing example is good. As an 
example consider the reference case of acceleration sensitivities calculated using 
the central difference method with updated modes shown in figure 5.50. For all 
derivative approximations, convergence is achieved with 32 or less modes. Modal 
convergence for acceleration sensitivities is equally good using the simple forward 
difference method with fixed modes as shown in figure 5.51. Figure 5.52 shows 
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the convergence of acceleration sensitivities calculated using the semi-analytical 
method with RWL vectors instead of vibration modes. Convergence is also good 
although slightly poorer than when modes are used. For example, approximately 
40 RWL vectors are required for a converged value of du t i p /dt J® compared with 
approximately 32 vibration mode shapes. 
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Figure 5.50. Modal convergence of tip acceleration sensitivities for the delta wing 
calculated using the central difference method. 


The modal convergence of stress derivatives, however, depends dramatically 
on whether fixed or updated modes are used in the calculation. The reference case 
with the central difference operator uses updated modes and as shown in figure 5.53 
the modal convergence for all stress sensitivities is very good. Also the convergence 
of the stress sensitivities using the forward difference operator with updated modes 
as shown in figure 5.54 is very good with 24 or less modes a yielding a converged 
solution. However, when the forward difference operator with fixed modes is used 
the modal convergence of the stress sensitivities is very poor as shown in figure 5.55. 
For one derivative approximation, da \\ the convergence is fairly good with 
approximately 24 modes yielding a converged solution. Especially poor convergence 
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Figure 5.52. 
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Modal convergence of tip acceleration sensitivities for the delta wing 
calculated using the forward difference method with fixed modes. 
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Modal convergence of tip acceleration sensitivities for the delta wing 
calculated using the semi-analytical method with RWL vectors. 
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is observed for dcr™ /dt\ 6 (derivative of stress in the lamina with respect to its own 
thickness) where approximately 100 modes are required for convergence. 
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Figure 5.53. 
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da«/dt ® (2.65) 
dt®/dt® (1.57) 


Modal convergence of selected stress sensitivities for the delta wing 
calculated using the central difference method. 


Using the semi- analytical method with fixed modes doesn’t improve the modal 
convergence of the stress sensitivities. Figure 5.56 shows the modal convergence of 
the same stress sensitivities as in the previous figures but calculated using the semi- 
analytical method with RWL vectors. The convergence behavior for each derivative 
approximation here is very similar to that for the forward difference method with 
fixed modes. 

However, when the basis vectors are assumed to vary with the design variables 
and the modified modal method (see section 4.2.2) is used to approximate d$/dx, 
the results are significantly improved. Figure 5.57 shows the modal convergence 
of the same stress derivative approximations as shown on previous figures. Here, 
the convergence is good with only around 24 modes required for convergence of the 
stress sensitivities. 
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Figure 5.54. 
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Modal convergence of selected stress sensitivities for the delta wing 
calculated using the forward difference method with updated modes. 
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Figure 5.55. 
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Modal convergence of selected stress sensitivities for the delta wing 
calculated using the forward difference method with fixed modes. 
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Figure 5.56. Modal convergence of selected stress sensitivities for the delta wing 
calculated using the semi- analytical method with RWL vectors. 
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Figure 5.57. Modal convergence of selected stress sensitivities for the delta wing 
calculated using the semi- analytical method with d^f/dx approxi- 
mated using the modified modal method. 
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It was mentioned in Chapter 4 that the modal method for approximating d$/ dx 
produces no improvement in the values of transient response sensitivities. This 
implies that including the modes in the modified modal method (see eq. (4.10) may 
also not significantly improve the transient response sensitivities. This implication 
was tested by studying the modal convergence of the stress sensitivities using the 
modified modal method but approximating d&/dx with only the pseudo-static term 
in equation (4.10). These results are shown in figure 5.58. Comparing this figure 
with figure 5.57 it can be seen that for more than 8 modes the results are nearly 
identical. It appears that a cheap, effective approximation to d$jdx in the semi- 
analytical formulation can be obtained using only the pseudo-static term from the 
modified modal method. 
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Figure 5.58. Modal convergence of selected stress sensitivities for the delta wing 
calculated using the semi- analytical method with d&/dx approxi- 
mated using only the pseudo-static solution. 


For the five-span beam example, the convergence of the stresses was substan- 
tially improved by including the static solution via either the mode acceleration 
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method, the static mode method, or the RWL method. The RWL method is at- 
tractive because it is cheaper to calculate n r RWL vectors than n r vibration mode 
shapes. However, incorporating the modified modal method in the sensitivity cal- 
culations with RWL vectors would seem to be impossible because it is derived to 
calculate the derivatives of vibration eigenvectors (see eq. (4.10)). Regardless, it 
seems like a worthwhile numerical experiment to try using RWL vectors along with 
the pseudo-static correction term from the modified modal method in the variable 
mode, semi-analytical formulation. One legitimate argument for doing this is the 
well-known observation that the basis spanned by the RWL vectors is an excel- 
lent approximation to the basis spanned by the eigenvectors. The results of this 
experiment for the modal convergence of the stresses in the delta wing are shown 
in figure 5.59. The convergence here is quite good also. For small numbers of 
modes the convergence is a bit erratic but in all cases the results are good for more 
than 32 modes. The benefit of combining the RWL vectors with the pseudo-static 
approximation to d§/dx is that the RWL vectors add the often important static 
displacement component to the basis while the pseudo-static term adds components 
reflecting the change in the design variable to the basis. 

As mentioned above, the benefit of the mode acceleration method is that it also 
includes this pseudo-static term. The semi-analytical, mode acceleration, sensitivity 
method described in Chapter 4 was also applied to this delta wing example. Again, 
the modal convergence of the stress sensitivities shown in the previous figures is 
considered. Figure 5.60 shows the excellent convergence of the stress sensitivities. 
Clearly, the mode acceleration method provides the same improvement in stress 
sensitivities as the semi-analytical method with a modified modal approximation to 
d$/dx. 
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Figure 5.59. Modal convergence of selected stress sensitivities for the delta wing 
calculated using the semi-analytical method with RWL vectors and 
the pseudo-static approximation to d$/dx. 
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Figure 5.60. Modal convergence of selected stress sensitivities for the delta wing 
calculated using the semi-analytical, mode acceleration method. 
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5.3. Stepped Cantilever Beam Example 


The third example considered is a cantilever beam with five, different rectan- 
gular cross sections along the length (see figure 5.61). This example is taken from 
ref. [38] where its minimum mass design under a static tip load was considered. 
The height and width of the beam cross section in each of the five sections are given 
in the table insert on figure 5.61 and represents an optimized design from ref. [38]. 
The beam is 200 inches long and, in the nominal case, each of the five sections has 
the same length. The material properties for the beam are also shown in fig. 5.61. 

In most of the analyses, the beam is modeled with three finite elements per 
section. The transverse displacement and rotation are the nodal unknowns resulting 
in a total of 30 degrees-of-freedom for this case. The effect of different numbers of 
elements per section on the the lowest 10 beam natural frequencies (with the beam 
clamped at the root) is shown in table 5.4. In the transient response analyses .5% 
of critical damping is included for each mode. 


Table 5.4. Lowest Frequencies for the Stepped Cantilever Beam 



Frequency, Hz 

Mode 

3-Elements 

4-Elements 

5-Elements 

6-Elements 

1 

22.67 

22.67 

22.67 

22.67 

2 

102.67 

102.66 

102.66 

102.65 

3 

249.72 

249.62 

249.80 

249.55 

4 

440.57 

440.04 

439.80 

439.67 

5 

652.50 

650.82 

650.04 

649.62 

6 

878.48 

874.48 

872.58 

871.54 

7 

1093.36 

1086.15 

1082.61 

1080.64 

8 

1296.61 

1285.63 

1279.94 

1276.72 

9 

1479.81 

1465.74 

1457.74 

1453.09 

10 

1641.44 

1625.75 

1615.41 

1609.19 


85 
















Section 

1 

2 

3 

4 

5 

Thickness, in. 

23.5 

22.0 

20.0 

18.0 

16.5 

Width, in. 

1.20 

1.10 

1.00 

.90 

.85 


Material Properties 
E = 30 x 10 6 lb/in. 2 

p = .3 lb/in. 3 
Length = 200 in. 

an applied rotational acceleration at 

5.3.1. Loading 

The loading for this stepped beam example is significantly different than for 
the first two examples. First, the load results from prescribing the acceleration 
at the beam root rather than by applying a force to the beam. Second, the time 
history as shown in figure 5.61 is more complicated than the simple step and ramp 
histories in the previous examples. The objective of this particular loading condition 
is to simulate the rotation of an appendage attached to a relatively large mass 
(e.g., robotic arm). The acceleration history in figure 5.61 rotates the root of 
the beam through 10 degrees in .18 seconds. After .18 seconds the beam root is 
motionless while other points in the beam are undergoing dynamic motion. Beam 
displacements, velocities, and accelerations in the following sections are with respect 
to the rotating coordinate system. 


Load history 
21.551 

0 0 

-21.55- 

.18 
t, secs 

Figure 5.61. Stepped cantilever beam with 
the root. 


rad 

sec 2 



This type of applied acceleration can be handled as an equivalent external force 
given as 

p = -Mr g(t) (5.1) 

where r is a vector describing the linear rigid body rotation of the beam about its 
root and g(t) is the prescribed acceleration history given in figure 5.61. It should 
be noted that the applied force in this case depends on the system mass matrix, 
and this must be considered in the sensitivity calculations. 

5.3.2. Stepped Beam Dynamic Response 

The transient behavior of the beam is strongly affected by the period of the 
loading. From table 5.4, the period of the lowest vibration mode is .044 secs while 
the period of the square-wave loading is .18 secs. From figure 5.62 it can be seen 
that in the time history of the beam tip displacement, this first mode predominates 
and almost exactly four cycles occur during the period of the loading. After the load 
is removed, the displacement response at the tip is relatively small. The bending 
stress at the root has a time history similar to that of the tip displacement as can 
be seen in figure 5.63 but with slightly more participation from higher frequency 
modes. As expected, the acceleration time history for the tip as shown in figure 5.64 
is considerably more jagged indicating the participation of many higher frequency 
modes. This behavior is largely due to the abrupt changes in loading in the square- 
wave input. Significant accelerations exist at the tip after the loading is removed. 
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Time, secs 

Figure 5.63. Time history of root stress for the stepped cantilever beam. 
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Figure 5.64. Time history of tip acceleration for the stepped cantilever beam. 

5. 3. 2.1. Modal convergence 

The first convergence study considered the effect of number of finite elements 
per section on the convergence of the critical point displacements, velocities, and 
accelerations at the beam tip and stresses at the root. For all these quantities the 
convergence is excellent. For example, the peak acceleration changes by less than 
one percent when the number of finite elements per section is varied from 3 to 8. 
Then, for the beam modeled with three elements per section, the effect of number 
of modes used in the analysis was considered. Generally the convergence was better 
than expected. Figure 5.65 shows the modal convergence of the tip acceleration at 
two different critical time points calculated using the mode displacement method. 
The values are essentially converged with five modes. 

The modal convergence of the stress at the beam root is also rapid. Figure 5.66 
shows the convergence of the root stress at two critical points calculated using the 
mode displacement method. No more than five modes are required for convergence. 
It was mentioned above that there is a strong static component in the beam response 
during the period while the load is applied. Usually this requires the use of the mode 
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Figure 5.65. Modal convergence of critical point tip accelerations for the stepped 
cantilever beam (mode displacement method). 


acceleration method or RWL vectors for acceptable convergence of the stresses. 
Evidently the lowest vibration mode is close enough to the static displacement shape 
for this cantilever beam that the mode displacement method gives good values for 
the stresses. 


5. 3. 2. 2. Using RWL vectors in the analysis 


In the stepped beam and delta wing examples the convergence with RWL 
vectors in analysis and sensitivity calculations was generally as good or better than 
with vibration modes. The modal convergence in the stepped cantilever beam 
example when RWL vectors are used is very good also as seen in figure 5.67 for 
accelerations. 

As can be seen in figure 5.67, the largest number of RWL vectors used in 
the analysis is 20. In the convergence studies considering vibration modes (e.g., 
fig. 5.65) the full set of 30 modes was used. A complete set of RWL vectors 
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Number of modes 

Figure 5.66. Modal convergence of critical point root bending stresses for the 
stepped cantilever beam (mode displacement method). 



Number of modes 

Figure 5.67. Modal convergence of critical point tip accelerations for the stepped 
cantilever beam (RWL vectors). 
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could not be generated for this example because of ill- conditioning inherent in the 
numerical process (eqs. (2.26), (2.27), (2.28), (2.29)). As additional RWL vectors 
are generated, roundoff errors cause the vectors to become less and less orthogonal. 
Eventually, the vectors become linearly dependent resulting in a singular reduced 
system. In most practical applications of this RWL method, this singularity problem 
would not occur because the number of RWL vectors generated would be much 
smaller than the total number of degrees of freedom. 

5.3.3. Sensitivities of Stepped Beam Dynamic Response 
5. 3. 3.1. Design variables 

Two different classes of design variables are considered in this example. The 
first class is the set of beam thicknesses in each of the five sections. They are 
denoted hi where i is the section number from figure 5.61. These are similar to 
thickness design variables considered in the five-span beam and delta wing examples. 
Sensitivity results are presented in the next sections considering hi and h$ from this 
set. 


The second class of design variables is the set of lengths of the five sections in 
the beam. The beam length is fixed at 200 inches; thus only four design variables 
determine the lengths of the five sections. The four design variables are denoted li 
where li is the distance from the beam root to the end of the ith section. Sensitivity 
results are presented in the next sections considering and I 4 from this set. In 
the structural optimization field this type of design variable is often referred to 
as a “shape” design variable and is studied separately from member thickness- 
type design variables. A recent study (ref. [39]) considered the calculation of 
static response sensitivities with respect to shape design variables using the semi- 
analytical method. It was found that numerical difficulties in the semi-analytical 
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method resulted in very large errors in sensitivities. This difficulty will be addressed 
in the next sections for the transient case. 

5. 3. 3. 2. Effect of finite difference step size 

It will be shown in this section that practically, the selection of finite difference 
step size is not a concern for this stepped beam example. A series of studies was 
performed to consider the effect of step size on both thickness and length sensitivities 
calculated using finite difference and semi- analytical methods. The finite element 
model with three elements per section was used and all 30 modes were included. 
Figure 5.68 compares approximate derivatives of tip displacement with respect to 
section thicknesses calculated using overall forward and central difference methods. 
A key point to be made is that both methods give excellent results for approximately 
an 8 decade step size range. For the large step size of .1 in., the central difference 
operator generally gives better results than the forward difference operator as would 
be expected. The results are nearly as good for sensitivities of the root stress with 
respect to the section thicknesses as shown in figure 5.69. Compared to figure 
5.68 there is slightly more error for the smallest and largest step sizes but the 
sensitivities are still accurate over a very broad range of step sizes. If sensitivities 
of stresses with respect to the length design variables are considered, the results 
are also very good. Figure 5.70 shows sensitivities calculated using the forward 
and central difference methods. Again there is a broad range of step sizes that 
provide accurate sensitivities. For the smaller step sizes the results are generally 
less accurate than in figure 5.69 but for the .01 step sizes they are more accurate. 

It is mentioned above that severe numerical difficulties were uncovered in ref. 
[39] when sensitivities of static response were calculated with respect to shape design 
variables. The result of this numerical ill-conditioning could be seen by calculating 
the sensitivities with different finite difference step sizes used for approximating 
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Figure 5.68. Effect of finite difference step size on the accuracy of tip displacement 
derivatives with respect to thickness design variables for the stepped 
cantilever beam (forward and central difference operators). 
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Figure 5.69. Effect of finite difference step size on the accuracy of root stress 
derivatives with respect to thickness design variables for the stepped 
cantilever beam (forward and central difference operators). 
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Figure 5.70. Effect of finite difference step size on the accuracy of approximate 
root stress derivatives with respect to length design variables for the 
stepped cantilever beam (forward and central difference operators). 


the derivatives of the stiffness matrix. For very small step sizes the error in the 
sensitivities is due to roundoff. For the larger step sizes, however, the errors in 
sensitivities were much larger than those due to truncation of the finite difference 
operator and were found to be caused by basic ill-conditioning in solving the semi- 
analytical equations. 


This same phenomenon occurs when sensitivities are calculated using a semi- 
analytical method for the transient case. Figure 5.71 shows approximate derivatives 
of root stress with respect to the length design variables calculated using the for- 
ward difference and semi- analytical methods. Again, all 30 modes are used in the 
analyses. For the smaller step sizes, the accuracy is significantly better for the 
semi-analytical method than for the overall forward difference method. For the 
10~ 2 step size, however, the results from the forward difference method are excel- 
lent while several of the sensitivities calculated using the semi-analytical method 
exhibit extremely large errors. This result is completely consistent with that in 
reference [39]. Although in this example there is a large range of step sizes where 
accurate sensitivities can be obtained, in general this would not be the case. Es- 
pecially as the problem becomes larger it is desirable to use larger step sizes in a 
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semi-analytical method but this is severely restricted for shape design variables by 
the type of error shown in figure 5.71. 
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Figure 5.71. Effect of finite difference step size on the accuracy of approximate 
root stress derivatives with respect to length design variables for 
the stepped cantilever beam (overall forward difference and semi- 
analytical methods). 


5. 3. 3. 3. Modal convergence of sensitivities 

Most of the sensitivities exhibit the same, good modal convergence as the re- 
sponse quantities. For example, the modal convergence of approximate derivatives 
of tip displacement with respect to hi and hs at different critical points are shown 
in figure 5.72. The sensitivities were calculated using the central difference method 
with updated modes and, as can be seen, the convergence is excellent. The conver- 
gence of tip acceleration derivative approximations is not as good as the displace- 
ment derivative approximations but is still acceptable as seen in figure 5.73. Again, 
these sensitivities are with respect to hi and /15 and calculated using the central 
difference method with updated modes. 

Convergence is also good when sensitivities with respect to the length design 
variables are considered. Figure 5.74 shows the modal convergence of approximate 
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Figure 5.72. Modal convergence of approximate derivatives of tip displacement 
with respect to thickness design variables for the stepped cantilever 
beam (mode displacement method, central difference operator). 
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Figure 5.73. Modal convergence of approximate derivatives of tip acceleration with 
respect to thickness design variables for the stepped cantilever beam 
(mode displacement method, central difference operator). 
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derivatives of acceleration with respect to l\ and U calculated using the central 
difference method. A step size of 10 -5 was used to avoid the problem shown in figure 
5.71. As can be seen in figure 5.74 convergence is achieved with approximately 10 
modes. 



Number of modes 

Figure 5.74. Modal convergence of approximate derivatives of tip acceleration with 
respect to length design variables for the stepped cantilever beam 
(mode displacement method, central difference operator). 


The modal convergence of stress sensitivities is similar to that for the delta wing 
example. When updated modes are used with an overall finite difference method, 
the convergence is excellent. An example of this is shown in figure 5.75 where 
approximate derivatives of the root stress with respect to hi and /15 calculated using 
the forward difference operator are shown. However, if fixed modes are used in a 
finite difference procedure, the modal convergence is much worse. Figure 5.76 shows 
an example of this for the same sensitivities as in figure 5.75. And if sensitivities 
of the root stress with respect to the length design variables (l lf I 4 ) are considered, 
the modal convergence is very poor. An example of this poor convergence is shown 
in figure 5.77. The convergence is similarly bad if the fixed mode semi-analytical 
method is used instead of a finite difference method. Figure 5.78 shows the poor 
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modal convergence of the same sensitivities as figure 5.77 calculated using the fixed 
mode, semi-analytical method. 



d °root /dh 1 (°- 02 > 
do root /dh 1 (°- 11 > 
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Number of modes 

Figure 5.75. Modal convergence of approximate derivatives of root stress with re- 
spect to thickness design variables for the stepped cantilever beam 
(mode displacement method, forward difference operator, updated 
modes). 


For the delta wing example, remedies for the poor convergence of stress sensitiv- 
ities in the semi-analytical method were based on approximating the mode shape 
derivatives, d&/dx. These semi-analytical methods including approximations for 
d#/dx were also applied to this stepped beam example. First d&/dx was approx- 
imated using the modified modal method. The modal convergence of the stress 
sensitivities is now excellent as can be seen in figure 5.79. The degree of improve- 
ment can best be appreciated by comparing figures 5.78 and 5.79 and noting that 
the range of the ordinate in figure 5.78 is much broader than in figure 5.79. Using 
only the first, pseudo-static term from the modified modal method as an approxi- 
mation to d&/dx was also tried. As can be seen in figure 5.80 the convergence is 
adequate though not quite as good as when the complete modified modal method 
is used. 
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Modal convergence of approximate derivatives of root stress with re- 
spect to thickness design variables for the stepped cantilever beam 
example (mode displacement method, forward difference operator, 
fixed modes). 
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Modal convergence of approximate derivatives of root stress with re- 
spect to length design variables for the stepped cantilever beam (mode 
displacement method, forward difference operator, fixed modes). 
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Figure 5.78. Modal convergence of approximate derivatives of root stress with re- 
spect to length design variables for the stepped cantilever beam (mode 
displacement method, semi-analytical method). 



Number of modes 

Figure 5.79. Modal convergence of approximate derivatives of root stress with re- 
spect to length design variables for the stepped cantilever beam (mode 
displacement method, semi-analytical, modified modal method). 
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Just as in the delta wing example, a case was also considered where RWL 
vectors were used instead of vibration modes but their derivatives were computed 
using the modified modal method (version with pseudo-static term plus modes). 
Again, somewhat surprisingly, the modal convergence of the stress sensitivities is 
good as seen in figure 5.81. 



Number of modes 

Figure 5.80. Modal convergence of approximate derivatives of root stress with 
respect to length design variables for the stepped cantilever beam 
(mode displacement method, semi-analytical, one-term modified mo- 
dal method). 


The semi- analytical, mode acceleration method was also tried as a remedy for 
the poor convergence of the stress sensitivities. Again, the very poor convergence 
is eliminated as can be seen in figure 5.82. 
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Number of modes 

Figure 5.81. Modal convergence of approximate derivatives of root stress with re- 
spect to length design variables for the stepped cantilever beam (RWL 
vectors, semi-analytical, modified modal method). 



Number of modes 

Figure 5.82. Modal convergence of approximate derivatives of root stress with re- 
spect to length design variables for the stepped cantilever beam (mode 
acceleration, semi-analytical method). 


103 



5.4. Summary 


A number of different methods for calculating sensitivities of transient response 
quantities have been exercised on three example problems, a five-span beam, a com- 
posite aircraft wing, and a variable cross section beam. Two of the methods are 
overall finite difference methods where the analysis is repeated for perturbed de- 
signs. The other methods are termed semi-analytical methods because they involve 
direct, analytical differentiation of the equations of motion with finite difference 
approximations of the coefficient matrices. All of the methods use basis vectors to 
reduce the dimensionality of the problem. Accordingly, the convergence of both the 
transient response quantities and their sensitivities as a function of number of basis 
vectors was a key concern in this chapter. 

In the delta wing and stepped cantilever beam examples the convergence of the 
response quantities was consistently very good. However, this was not the case with 
the five-span beam. With the five- span beam under a concentrated end moment and 
ramp time history the convergence of displacements and velocities was adequate. 
However, the convergence of accelerations was poor. The convergence of stress 
resultants for this example depended on how they were calculated. When the mode 
displacement method was used, the convergence was quite poor. However, when the 
mode acceleration method, the Ritz-Wilson-Lanczos vector method, or the static 
mode method was used, the convergence was good. In cases where convergence was 
poor for the five-span beam, the addition of modal or discrete damping improved 
the convergence somewhat. However, it did not eliminate the convergence problems. 

The modal convergence of the sensitivities in the three examples is consistent 
with the convergence of the response quantities themselves. For the delta wing and 
stepped cantilever beam examples the convergence of sensitivities was generally 
good. For the five-span beam example the convergence of displacement sensitivities 
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was adequate but the convergence of velocities, accelerations, and stress resultants 
was generally poor. This poor convergence was observed for all the sensitivity 
calculation methods. Furthermore, it appears to be associated with the structure 
and loading because no improvement was observed as the number of finite elements 
per span was increased. 

In certain cases poor convergence of sensitivities was also observed for the 
delta wing and stepped cantilever beam examples. When sensitivities of stresses 
were calculated using the fixed mode, overall finite difference methods or the semi- 
analytical method assuming fixed modes, the convergence was very poor. In large 
problems, however, updating the vibration modes in the overall finite difference 
methods or rigorously calculating derivatives of the mode shapes is very expensive. 
The mode acceleration version of the semi-analytical method and the semi-analytical 
method with mode shapes approximated using the modified modal method were 
devised to alleviate this poor convergence with lower computational cost. When 
both methods are applied to delta wing and stepped cantilever beam examples the 
modal convergence of sensitivities is excellent. 

All of the sensitivity calculation methods considered herein rely on finite dif- 
ference operators. Thus step size selection is an important concern. The system 
stiffness and mass matrices are linear functions of many of the design variables in 
the three example problems. This allowed large step sizes to be used in the semi- 
analytical methods to minimize the roundoff errors and produce accurate derivatives 
of the stiffness and mass matrices. Also there is less opportunity for roundoff error 
in calculating finite difference derivatives of just the coefficient matrices compared 
with finite difference derivatives of the overall response quantities. For these rea- 
sons, the semi-analytical methods were consistently less sensitive to finite difference 
step size than the overall finite difference methods. 
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Chapter 6 


Computational Costs 

A consideration of the computational costs is essential for evaluating any nu- 
merical method. This is especially difficult in large-scale, finite element-based pro- 
cedures because there is often considerable “overhead” required in the practical 
implementation of a given numerical method. For example, most finite element 
codes require only a small portion of a system matrix to be resident in central 
memory during factorization at any given time. The other portions of the matrix 
are read from disk and the factored portions are written to disk as required. A 
similar situation can exist on virtual memory machines where the disk operations 
are transparent to the implementor. In these cases, the computer resources required 
are very implementation dependent. 

An approach that is common in the formal study of numerical methods is to 
evaluate the computational cost by counting the number of floating point oper- 
ations. There are some pitfalls to this approach. In some cases, even for large 
problems, because of the required overhead it is impossible to achieve a practi- 
cal implementation that will execute as fast as the predictions from the operation 
count. In other cases, especially on vector machines, it is possible for a method 
with a higher operation count to be faster than a method with a lower operation 
count. 

Nevertheless, this approach will used here, primarily to indicate the major 
trends in the costs of the methods but not to make fine distinctions between them. 
Following common practice, a floating point operation (or “flop”) is defined as the 
combination of a floating point multiply, add, and associated array indexing. In the 
remainder of this chapter a floating point operation is often referred to simply as 
an operation. 
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6.1. Costs of Basic Matrix Manipulations 


Multiplication of full matrices occur in several places in the transient response 
and sensitivity methods. The approximate number of floating point operations 
required to multiply a full l x m matrix and a m x n matrix is given as 


Cfmul = l mn (6.1) 

Solution of the reduced eigenproblem (eq. (2.23)) is important in solving the 
system eigenproblem using subspace iteration (eq. (2.3)) and in uncoupling the 
reduced system when basis vectors other than the eigenvectors are used. In both 
cases it is necessary to solve a full, generalized eigenproblem for all n r eigenvalues 
and eigenvectors. Since eigenvalue solution techniques are inherently iterative, the 
number of operations required for a converged solution can only be estimated. Ref- 
erence [15] estimates the number of operations for complete solution of a generalized 
eigenproblem using the Jacobi method as 


C reig = 18n* + 36n" (6.2) 

Other techniques for solving this eigenproblem may have a significantly different 
cost. However, it will be shown that the cost of this eigensolution is small relative 
to other tasks in the sensitivity calculations. 
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0.2. Costs of System Matrix Manipulations 


For the purpose of considering the computational costs of operations on system 
matrices (e.g., K, M), these matrices are considered to be stored in a banded form. 
In a banded form, only matrix elements located near the diagonal are stored; the 
matrix elements outside this “bandwidth” of the matrix are zero and are not stored 
or considered in operations. Finite element problems yielding a stiffness matrix with 
a constant bandwidth are rarely encountered in practice so most finite element 
codes use more sophisticated and efficient schemes for storing system matrices. 
However, having a single, easily-understood number to characterize the sparsity 
of a system matrix (the bandwidth) is convenient in approximating computational 
costs. Although few finite element problems have precisely a constant bandwidth, 
this assumption is accurate enough in many cases to get reasonable estimates for 
relative number of operations in a numerical procedure. 

From ref. [40] the cost in number of operations of factoring a banded system 
matrix of order n g is given as 


= i/?(/3 + 3)n, - y - /3 2 - f/» (6.3) 

where (3 is half the bandwidth (excluding the diagonal) Also from ref. [40], the cost 
of a single solution of a banded system, given the factored matrix, is given as 


Cbtol = 2(/3 + l)n 9 — /?(/? + 1) 


(6.4) 


The cost of multiplying a banded system matrix and a single vector is given as 


Cbmul — (2/3 + 1)% 


( 6 . 6 ) 
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6.3. Cost of Basis Reduction 


The process of reducing the degrees of freedom from n g to n r requires the 
matrix triple product operations shown in equations (2.6), (2.7), and (2.8). Since 
this process is used in all the sensitivity methods, the cost will be considered sep- 
arately here. In performing this operation, n r system vectors are multiplied by a 
system matrix. Then n r (n r + l)/2 inner products (for a symmetric system matrix) 
between system vectors are performed. The total number of required floating point 
operations is 

Cred = n r Cbmul + Tl g Tl r (n r + l)/2 (6.6) 

6.4. Cost of System Eigensolution 

The cost of solving the generalized vibration eigenvalue problem is even more 
difficult to estimate than the cost of solving the reduced eigenproblem. The nu- 
merical techniques vary widely among different analysis codes. Furthermore, a 
technique used for one problem might be totally inappropriate when applied to 
a different problem. Nevertheless, some assumptions will be made here that will 
hopefully lead to a reasonable estimate of computational costs for a fairly broad 
class of problems. 

First, it will be assumed that the eigenvalue problem will be solved using a 
subspace iteration technique with shifts (see for example ref. [15]). In recent years, 
software based on this approach has become common. Also, the eigensolver, E4, 
in the EAL software used in this study (ref. [23]) is based on this approach. It 
is also necessary to make assumptions about the number of vectors used in the 
subspace and the number of iterations at a given shift point required to converge 
some subset of these vectors to eigenvectors. The following numbers were used 
for these quantities with the realization that they may be optimum for only a few 
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problems. Also, it is assumed that the eigensolution is being performed for a slightly 
perturbed model and the eigenvectors from the initial model are available as the 
initial subspace. At each shift point = 16 vectors are included in the subspace. 
After Tin = 2 iterations, n caa = 8 of these vectors have converged to eigenvectors. 

The number of shifts or number of factorizations required will be approximately 

n„Kift = -^ !L + 1 (6.7) 

n ct t 

At each iteration, the inverse power operation requires a matrix product between 
M and nt aa vectors followed by nt aa solutions of the system equations based on 
the current, factored K. The basis reduction operation for both K and M requires 
TTt aa {n taa + 1) system vector inner products. Next, the eigenproblem of reduced 
order n*,, must be solved. This cost is given in equation (6.2) with n r replaced 
by nt aa . Finally, the updated set of approximate, system eigenvectors must be 
formed as a linear combination of the current approximation. This requires n\ at n g 
operations. The approximate cost of solving the system eigenproblem for n r modes 
and frequencies can be written as 

C e ig = Tl a hiftCbfac'^~ 

Tl a hiftTlit(^Tlt aa Cbmul 4" TltasCbsol 4“ Tlt aa (^Tlt aa 4" l)n.g (6*6) 

4- 18n L + 36 ™L + n u, n 9 ) 

6.5. Cost of Generating RWL Vectors 

As has been demonstrated, RWL vectors are an attractive alternative to vi- 
bration mode shapes for basis reduction in transient response analysis. It has been 
mentioned previously that generation of the RWL vectors is considerably cheaper 
than vibration modes. An estimate of this cost in number of floating point opera- 
tions will be derived here. 
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First a factorization of the system K is required. The system equations are 
solved n r times based on the factored K. The generation of right-hand-side vectors 
requires n r — 1 matrix products between M and a vector. Another key step in the 
process is the Gram-Schmidt orthogonalization as indicated in equation (2.27). For 
all vectors, this requires n r — 1 multiplications of a vector by M and n r (n r — l)/2 
vector inner products. The scaling of each vector requires n r vector inner products 
and n r divisions of a system vector by a scalar. Writing the total number of floating 
point operations in expanded form yields 

Crwl = Cbfac Cb>oi + 2(n r — l)(7& m ,w H — ^ - n g + 2n r n 9 (6.9) 

6.6. Cost of Model Generation 

The generation of the finite element model requires processing of the input, 
forming elemental matrices, and forming global, system matrices. Most of the sen- 
sitivity calculation methods require generation of a single perturbed model for each 
design variable. The central difference method, however, requires the generation 
of two perturbed models. Thus to compare the central difference method with the 
other methods, an estimate of the model generation cost is required. This cost is 
difficult to calculate in general. For the purposes herein this cost is estimated em- 
pirically with EAL by observing the execution time for model generation relative 
to matrix multiplication for a number of models. From these experiments it was 
observed that the predominant element type in the model substantially affects the 
cost. That is, forming the element matrices in a model composed of three dimen- 
sional solid elements is much more costly than in a model composed of rod elements. 
The estimate for model generation cost used here, 

Cmodel = 100 /3n g (6.10) 

roughly approximates the cost for a model with 2-D, plate-type elements in EAL 
but would be significantly in error for predominantly 1-D or 3-D models. 
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6.7. Cost of Integration of Reduced System 


The basic operation for integrating the reduced system is shown in equation 
(2.30). The two matrix multiplications shown in equation (2.30) are performed at 
every time step. If equations (2.5) are coupled, the W ij and N^- are full and the 
explicit matrix multiplication must be performed. In this case, the total number of 
floating point operations for integration of the system is given as 

Cinte = 8n*n t (6.11) 

where n* is the number of time steps in the analysis. If equations (2.5) are not 
coupled, the W ij and N ij are diagonal and this fact can be exploited to substantially 
reduce the cost of integrating the system. The number of floating point operations 
in this case is given as 

Cinte 8n r n t (6.12) 

When the number of equations in the reduced system, n r , is large, the difference 
between Cinte in equations (6.11) and (6.12) is very large. For the comparisons of 
sensitivity methods in this chapter, equation (6.12) is used to estimate the integra- 
tion cost. When vectors other than vibration modes are used or vibration modes for 
an initial model are used with a perturbed model, the equations are first uncoupled 
by solving the reduced order eigenproblem. 

6.8. Cost of Back Transformation for Physical Response Quantities 

After the reduced equations have been solved, it is necessary to recover the 
physical displacements, velocities, accelerations, and stresses (or stress resultants) of 
interest. Usually the quantities of interest are only a subset of all possible quantities 
available from the finite element model. In the critical point constraint formulation 
described in Chapter 3 it is necessary to recover the physical response quantities 
only at the critical times. That is, the back transformation is performed at only 5 


112 



to 20 critical points rather than at thousands of time steps. The cost of the basic 
back transformation operation is 


C back — ft’p Tlr (6.13) 

where n p is the number of physical quantities being recovered from the modal values 
and n c is the number of critical points. The costs of back transformation in the 
specific sensitivity methods will be expressed as a multiple of this basic cost. 

0.9. Cost of Sensitivity Calculation Methods 

Because all of the sensitivity calculation methods require the dynamic analysis 
of the initial model, this component of the cost can be neglected in comparing the 
different methods. In addition, in all the methods the basic operations are repeated 
for each design variable so the costs estimated below are per design variable. Also, 
to simplify the cost analysis, the models are assumed to be undamped so that 
any operations dealing with modal damping or system damping matrices are not 
included. 

6.9.1. Finite Difference Methods 


Both forward and central difference methods for calculating sensitivities were 
considered in Chapter 4. In the central difference method, the basic operations of 
the forward difference method are performed twice, so the cost is approximately 
twice that of the forward difference method. Costs will be derived here for the 
forward difference method. In both finite difference methods, the basis vectors can 
be the same as for the original model (fixed) or recalculated for the perturbed model 
(updated). The cost with updated modes presented here is based on using natural 
vibration modes. An alternative of using RWL vectors is considered separately. 
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The first step in the forward difference method is evaluation of the perturbed 
model. If the modes are being updated, the eigenproblem is solved. Otherwise, the 
original modes are used to reduce the basis and the reduced order eigenproblem 
solved to uncouple the transient equations. The uncoupled equations for the per- 
turbed system are then integrated and the n p physical quantities calculated at the 
n c critical time points. The cost of the actual difference operation is very small and 
is therefore neglected. For the fixed mode case, the total cost is 

C fdfix — C model "b ^C re d "I" Creig "b d" C% n te "b Chuck (6.14) 

For the updated mode case, the total cost is 

C fdupd — Cmodel "b C e ig "b C inte + Cl, ack (6.15) 


6.9.2. Semi-Analytical Method With Fixed Modes 

The semi-analytical method begins by evaluating the perturbed model. Then 
dM/dx and dK/dx are formed using a forward difference operator. Each derivative 
requires about (3n g operations. Then the basis reduction operation is applied to both 
derivative matrices. Formation of the right-hand-side pseudo load (see equation 
(4.6)) is a fairly costly operation and the two matrix products dM/dx q and dfC/dx q 
require about n 2 r nt operations each. Finally the uncoupled equations are integrated 
and the physical sensitivities recovered. For the purposes of cost estimation a single 
quantity n p is used as the total number of required physical sensitivities. In the semi- 
analytical methods, however, the specific procedure for recovering the sensitivities 
depends on whether the quantity is a displacement, velocity, acceleration, or stress 
sensitivity. In estimating the costs of this back transformation operation these 
differences are ignored. One justification for this approach is that the cost of back 
transformation is usually small relative to other costs in the sensitivity calculation. 
In this fixed mode, semi-analytical method, approximately the same number of 
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operations are required for the recovery of physical sensitivities as in the finite 
difference methods. The total number of floating point operations can be written 
as 

C safix — Cmodel 4" 4" 2 Cred 4" 2,7l r TLi Cinte + Cback (6.16) 

6.9.3. Semi- Analytical Method With Approximate d&/dx 

Just as in the fixed mode semi- analytical method, evaluating the perturbed 
model and forming dM/dx and dK/dx is the first step. The next step is using the 
modified modal method to approximate d$/dx. 

The procedure for the modified modal method is given in equations (4.9), 
(4.10), (4.11), and (4.12). The calculation of the n r pseudo-static contributions 
requires the formation of n r right-hand-side vectors and n r solutions of the system 
equations. The formation of the Ajk participation factors requires approximately 
n r system matrix additions plus the equivalent of a triple product basis reduction 
operation. Forming the linear combination of pseudo-static term and eigenvectors 
requires n r n g operations. The total cost in number of floating point operations for 
the modified modal method is 


Cmmod — 2n r Tlg -f- Tl r Cbsol 4" 4" Cred 4“ Tl r Tlg (6.17) 

Given d&/dx, the derivatives of the reduced system matrices can be formed. 
For both dM./dx and d&/dx two triple product, basis reductions plus n r vector 
inner products (for the d& T /dx M# term since M4> is already available) are re- 
quired as shown in equation (4.8). The right-hand-side formation and integration 
of the reduced sensitivity equations are identical to the fixed mode semi- analytical 
method. Because of the non-zero d&/dx, recovery of the physical sensitivities is 
more complicated than in the fixed mode case. Approximately twice the number of 
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operations are required in the back transformation since both $ and d&/dx terms 
must be considered as shown in equations (4.13). The total cost for the variable 
mode semi-analytical method is 


Csaupd — Cmodel 4“ ‘ifitlg -f- C rnmod "1" 4 C r ed 4" 2 Tl r Tlg -|- 2Tl r 7lf 4" Ci n te 4“ 2(7 back (6.18) 


6.9.4. Semi- Analytical. Mode Acceleration Method 

Since di\/dx and dq/dx are obtained from the fixed mode semi-analytical 
method, the operations in equation (6.16) (except Cback) are required in applying 
the mode acceleration method. The back transformation operations for displace- 
ment and stress sensitivities are more complicated as seen in equations (4.19) and 
(4.20). The cost of forming the coefficients in equation (4.19) is dominated by mul- 
tiplying a vector by a system matrix, adding n r + 1 system vectors, and solving the 
system equations for n r + 1 pseudo-static vectors. Again, the assumption is made 
that the model is undamped so the C and the dC/dx terms in equation (4.19) are 
zero. 


The back transformation procedure for displacement and stress sensitivities 
involves application of equations (4.19) and (4.20) for each quantity at the crit- 
ical times. Velocity and acceleration terms are calculated as in the fixed mode, 
semi-analytical method. Again, using only a single quantity for the number of 
back transformed quantities n p , the cost can only be roughly estimated as 4 Cback- 
The total cost for the semi-analytical, mode acceleration method can then then be 
written as 

C Mamacc —C model 4" 2/?Tlg -f- 2 Cred 4” 2tl r Tlf 

(6.19) 

4" Cbmul 4" (Wf 4- 1)^5 4“ (Wf 4* 1 'jCbsol 4- Cinte 4" 4 Cback 
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6.10. Analysis of Cost For Various Models 


With the above expressions for computational cost it is now possible to evaluate 
the use of the sensitivity calculation methods on various examples. The first three 
examples are those considered in Chapter 5. These three examples, however, are 
all rather small compared with the class of problems envisioned for the production 
use of the sensitivity methods. Accordingly, two other hypothetical problems with 
a larger number of degrees of freedom have been included. 

The key parameters from the five problems required for the cost analysis are 
shown in table 6.1. Several points should be made about these parameters. The 
two beam problems have a small number of degrees-of-freedom and a very small 
bandwidth and, as a result, a small cost for system matrix factorization. This is 
unusual in finite element analysis. Medium model A and large model B represent a 
typical medium size, linear dynamics problem and a rather large, ambitious problem 
respectively. Medium model A also is complicated by the fact that 100 vectors are 
assumed to be required in the transient analysis. In all five examples, a relatively 
large number of time steps are used in the transient analysis. 


Table 6.1. Parameters Governing Computational Costs 


Models 

n 9 

/3 

n r 

n t 

n p 

n c 

Five-span-beam 

32 

3 

18 

6000 

25 

10 

Delta Wing 

264 

30 

20 

30000 

13 

5 

Stepped Beam 

32 

3 

20 

30000 

4 

5 

Medium Model A 

3000 

100 

100 

10000 

50 

10 

Large Model B 

12000 

300 

30 

20000 

200 

10 
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6.10.1. Cost of Computational Subtasks 


Table 6.2 shows the number of floating point operations required for different 
computational tasks for the five example problems. Examining the costs for these 
subtasks gives some clues to the costs of different sensitivity calculation methods. 
For the first three examples the cost of system matrix factorization is low. For the 
two larger, hypothetical examples the factorization cost is much higher relative to 
that of other tasks. In the first three examples, the cost of integrating the reduced 
equations is substantial even though the equations are uncoupled. For the models 
A and B, the integration cost is one to two orders of magnitude less than the other 
subtask costs in table 6.2. Consistently, in all five examples, the cost of performing 
the triple product, basis reduction is high. For the three small problems this cost is 
significantly higher than the factorization cost. For medium model A this cost is also 
much higher than the factorization cost but this is primarily due to the requirement 
of 100 vectors in the reduced system. Even in model B, however, the basis reduction 
cost is only a little less than half the factorization cost. One conclusion of this is 
that the number of vectors in the reduced system substantially affects the cost of 
the analysis even if the vectors are not updated for the current model. 

The use of RWL vectors in the transient and sensitivity analyses was considered 
in Chapter 5. Here, the cost of generating RWL vectors compared with vibration 
modes will be considered. Table 6.2 shows the cost of system matrix eigensolution, 
C e ig and RWL vector generation, C rw [ for the five example problems. In every 
case the generation of RWL vectors is cheaper than the eigensolution. In the beam 
examples C rw i is more than an order of magnitude less than C'ig- This results from 
the unusual situation in which the number of required eigenvectors is nearly the 
same as the total number of degrees of freedom. In this case, the solution of the 
reduced eigenproblem artificially raises the cost of the system eigensolution. The 
other three examples show C e i g to be 3 or 4 times C rw l- This is probably a much 
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more accurate estimate of the cost savings obtained by using RWL vectors instead 
of eigenvectors. 


Table 6.2. Number of Operations for Selected Computational Subtasks 


Models 

Cbfac 

Cred 

Ceig 

Crwl 

Cinte 

Five-span-beam 

2.7 x 10 2 

9.5 x 10 3 

6.6 X 10 5 

1.8 x 10 4 

8.6 x 10 5 

Delta Wing 

1.2 x 10 5 

3.8 x 10 s 

4.8 x 10 6 

1.1 x 10 6 

4.8 x 10 6 

Stepped Beam 

2.7 x 10 2 

1.1 x 10 4 

6.6 x 10 5 

2.1 x 10 4 

4.8 x 10 6 

Medium Model A 

1.5 x 10 7 

7.5 x 10 7 

7.4 x 10 8 

2.1 x 10 8 

8.0 x 10 6 

Large Model B 

5.4 x 10 8 

2.2 x 10* 

4.0 x 10 9 

1.2 x 10 9 

4.8 x 10 6 


6.10.2. Comparison of Costs for Five Sensitivity Methods 


The primary objective of this chapter is to compare the costs in number of 
floating point operations of the sensitivity methods. This is summarized for five 
sensitivity methods, for the five examples in table 6.3. It is believed that these 
five sensitivity methods are all practical alternatives for large-order problems. This 
belief is substantiated by the fact that for all five examples the difference among 
the five costs is less than one order of magnitude. 


Table 6.3. Overall Operation Costs for Five Sensitivity Methods 


Models 

H2S59H 

M3R 

mm 

■3 

i 

H 

■3 

BH 

Ctamacc 

Five-span-beam 

1.0 

X 

10 6 

1.5 

X 

q m 

4.8 

X 

10 6 

4.8 

X 

To^ 

4.8 

X 

10 6 

Delta Wing 

6.5 

X 

10® 

1.0 

X 

yjf 

3.0 

X 

10 7 

3.2 

X 

10 7 

3.1 

X 

10 7 

Stepped Beam 

5.0 

X 

10 6 

5.5 

X 

\m 

2.9 

X 

10 7 

2.9 

X 

10 7 

2.9 

X 

10 7 

Medium Model A 

2.4 

X 

10 8 

7.8 

X 

10 8 

3.9 

X 

10 8 

6.8 

X 

10 8 

4.5 

X 

10 8 

Large Model B 

8.2 

X 

10 8 

4.4 

X 

10® 

8.5 

X 

10 8 

1.7 

X 

10® 

1.1 

X 

10® 


The forward difference method with fixed modes is consistently the cheapest 
method. However, this low computational cost must be weighted against the pitfalls 
of the method discussed in Chapter 5. The cost of a fixed mode central difference 
method which is approximately twice the forward difference cost would also be 
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quite competitive with the other methods and would lessen the sensitivity to finite 
difference step size. For the two larger problems, the forward difference method 
with updated modes is relatively expensive and an updated mode central difference 
method would be extremely expensive for larger problems. 

In the three smaller problems the semi- analytical methods require significantly 
more operations than the finite difference methods. This is primarily because the 
larger number of time steps makes the calculation of the right-hand-side pseudo- 
load relatively large. For the two larger problems, however, the fixed mode semi- 
analytical method is quite competitive with the forward difference method. For 
model A it is less than twice the cost of the finite difference method and for model 
B it is essentially the same. 

The number of basis vectors used is a key parameter in both the analysis 
and sensitivity calculations. Table 6.1 shows the number of modes used in the 
baseline cost analyses for the five examples. Here, the effect of number of modes 
on the overall sensitivity costs is considered. First, the delta wing example, which 
is representative of a typical small problem is considered. Shown in figure 6.1 is 
the cost for the five methods plotted as a function of number of vectors used. The 
number of modes ranges from 20 to 100. The values of the other parameters in the 
problem are in table 6.1. The key result from figure 6.1 is that the semi-analytical 
methods are much more costly than the finite difference methods for large numbers 
of modes. There are two reasons for this. First, because the problem is small, 
calculation of the vibration modes is relatively cheap. Second, because there are 
a large number of time steps, formation of the right-hand-side in the sensitivity 
equations for the semi-analytical methods is quite costly when the number of modes 
used is large. 

For the large model B example, the result of varying the number of modes is 
very different. For this example the cost of the five sensitivity methods plotted as a 
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function of number of modes is shown in figure 6.2. In this example, the calculation 
of the modes is a very costly operation. Accordingly, the forward difference method 
with updated modes is substantially more costly than the other methods for large 
numbers of modes. The fixed mode forward difference and semi-analytical methods 
show only moderate increases in cost as the number of modes is increased. 

It was mentioned above that a relatively large number of time steps are used 
in the five examples. The effect of number of time steps on the overall sensitivity 
calculation costs is considered here. The delta wing example is considered as rep- 
resentative of a small problem and the large model B example is representative of 
a large problem. For the delta wing the computational costs for the five sensitivity 
methods are plotted sis a function of number of time steps in figure 6.3. The values 
of the sensitivity calculation costs here are similar to those in figure 6.1; the for- 
ward difference methods show only moderate cost increases for larger numbers of 
time steps and the semi-analytical methods show substantial cost increases. Again, 
the reason is that the right-hand-side formation in the semi-analytical methods is 
a substantial part of the total cost in small problems. 

The cost results from a large problem, the model B example, are very different, 
however. The costs as functions of number of time steps are plotted in figure 
6.4. Obviously, there is practically no change in cost for any of the methods as a 
function of number of time steps. For this large problem both the cost of integrating 
the uncoupled equations and the cost of forming the right-hand-side in the semi- 
analytical methods are two to three orders of magnitude less than the total cost. 
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Figure 6.4. Cost of the sensitivity calculation methods as a function of number of 
time steps for the model B example. 
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6.11. Summary 


The main objective of this chapter is a comparison of the computational costs 
in number of floating point operations of the sensitivity calculation methods. Five 
example problems were considered — the three example problems from Chapter 5 
which are all fairly small and two larger, hypothetical examples. 

Many of the results depend significantly on whether the problem is one of the 
three smaller examples or one of the two larger, hypothetical examples. In the three 
smaller examples, the cost of system matrix factorization is low while in the larger 
problems this cost is quite high. When the cost of factorization is high the system 
eigenproblem is especially costly. In the smaller problems, operations repeated 
for the reduced problem at each time step (such as integration of the uncoupled 
equations) are a significant percentage of the total sensitivity calculation cost. For 
large problems the relative cost of these operations is small. 

For all five examples the forward difference method with fixed modes was the 
cheapest. For the smaller problems the forward difference method with updated 
modes had a relatively low cost but for the larger problems the cost was quite high. 
For the larger problems the semi- analytical method with fixed modes and the semi- 
analytical, mode acceleration method have costs that are relatively competitive with 
the fixed mode forward difference method. In till cases, the semi-analytical method 
with approximate eigenvector derivatives was one of the more costly methods. 

It was shown in Chapter 5 that for two examples the accuracy of the stress sensi- 
tivities for small numbers of basis vectors was extremely poor. It was demonstrated 
that the semi- analytical, mode acceleration method was one means of dramati- 
cally improving this accuracy. From the results of this chapter, the semi-analytical, 
mode acceleration method is only slightly more costly than the fixed mode forward 
difference and semi- analytical methods. Given the unacceptable accuracy of these 
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fixed mode methods for two of the examples, the semi-analytical, mode acceleration 
method appears to be the method of choice. 
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Chapter 7 


Concluding Remarks 

Several methods have been developed and evaluated for calculating sensitivities 
of displacements, velocities, accelerations, and stresses in linear, structural, tran- 
sient response problems. Two of the methods are overall finite difference methods 
where the analysis is repeated for perturbed designs. The other methods are termed 
semi-analytical methods because they involve direct, analytical differentiation of the 
equations of motion with finite difference approximations of the coefficient matri- 
ces. The different sensitivity methods were evaluated by applying them to three 
example problems, a five-span, simply supported beam loaded with an end moment, 
an aircraft wing loaded with a distributed pressure, and a cantilever beam with a 
stepped cross section loaded with an applied root angular acceleration. 

An important issue in calculating transient response sensitivities for use in 
formal optimization procedures is how to define the constraints. Two common 
approaches are to integrate the response quantity over time or to pick the maximum 
(or minimum) value of the response quantity in time. Both of these approaches 
have drawbacks. An alternative critical point constraint approach was implemented 
which identifies the most important response points along the time history. A 
method for identifying these critical points was devised that, based on the three 
examples considered, appears to be very effective even for very jagged response 
histories. 

All of the analyses and sensitivity methods considered use approximation vec- 
tors to reduce the number of degrees of freedom in the analysis. Vibration mode 
shapes, Ritz- Wilson- Lanczos vectors, and static displacement shapes were used in 
the analysis and sensitivity calculations. The key question when an approximate, 
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reduced basis is used in an analysis is how many basis vectors are required for an 
accurate approximation to the finite element solution. It was generally found that 
if the accuracy of the response quantities was poor, the accuracy of the sensitivities 
was extremely poor. In a number of cases, however, even though the accuracy of the 
response quantities was adequate, the accuracy of sensitivities was poor. This will 
be discussed further below. In all cases considered here, the accuracy as a function 
of number of vectors for both the response quantities and sensitivities with Ritz- 
Wilson-Lanczos vectors was as good or better than with vibration modes. Since the 
generation of Ritz-Wilson-Lanczos vectors is cheaper than vibration modes, they 
appear to provide a more cost effective alternative to modes in many cases. 

A goal in considering sensitivity methods in this study is that they be suitable 
for very large-order finite element analysis. In these types of problems a complete 
vibration analysis for each perturbed model is impractical because of the high com- 
putational cost. To reduce this cost, one approach which was studied herein is to 
use the basis vectors from the initial model to approximate the response in the 
perturbed model. This often provides an effective solution. In two of the three 
examples problems considered, however, using the initial vectors in an overall finite 
difference method or assuming fixed modes in a semi-analytical method resulted in 
very poor modal convergence for stress sensitivities. Two methods were devised to 
improve this poor performance. 

The first method retains the derivatives of the basis vectors in the sensitivity 
equations but approximates these derivatives rather than using a very costly exact 
computation. One well-known method for approximating eigenvector derivatives, 
the modal method, was found to be completely ineffective because it adds no new 
information to the existing modal basis. Another technique, the modified modal 
method, adds a pseudo-static contribution to the eigenvectors in approximating the 
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eigenvector derivatives. This technique, along with the semi-analytical method, was 
found to be very effective at improving the poor accuracy of the stress sensitivities. 

A second method for improving the accuracy of the stress sensitivities as a 
function of number of modes is to use a mode acceleration version of the semi- 
analytical method. The key to the mode acceleration method in the transient 
analysis is that it supplements the modal basis with a static contribution calculated 
from the complete model. The key to the mode acceleration implementation of 
the semi-analytical sensitivity method is that it supplements the modal basis with 
pseudo-static sensitivity terms calculated from the complete model. This technique 
produced the same dramatic improvement in the accuracy of stress sensitivities as 
the semi- analytical, modified modal method. 

As mentioned above, computational cost was an overriding concern in consid- 
ering the sensitivity analysis methods. To estimate this cost, expressions for the 
number of floating point operations in each of the methods were derived. Although 
this approach doesn’t include important effects such as overhead operations or disk 
I/O that would be present in a practical implementation of these methods, it does 
a provide a mechanism for an approximate, coarse ranking of the methods by com- 
putational cost. The overall forward difference method with fixed basis vectors was 
found to be the cheapest method for all cases considered. This technique, however, 
suffers from the accuracy problems mentioned above. One approach to alleviat- 
ing these accuracy problems is to recalculate the modes for the perturbed model 
(updated modes) in the overall forward difference method. This forward difference 
method with updated modes was found to be very costly for large, models, however. 
The fixed mode, semi-analytical method is only slightly more costly than the over- 
all forward difference method with fixed modes, but suffers from the same accuracy 
problem as the fixed mode, overall forward difference method. Two techniques with 
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reasonable costs that alleviate the accuracy problem are the mode acceleration im- 
plementation of the semi-analytical method and the semi-analytical method with 
approximate mode shape derivatives. Of these two methods, the semi-analytical, 
mode acceleration method is slightly cheaper. 

Given the high accuracy of the semi- analytical, mode acceleration method for 
a relatively small number of modes and its reasonable computational cost, this 
appears to be the method of choice. In the three examples considered herein, this 
method consistently performed as well as the much more costly, updated mode, over- 
all finite difference methods. Furthermore, the insensitivity of this and the other 
semi-analytical methods makes this semi-analytical, mode acceleration method es- 
pecially attractive. 
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Appendix 


Computer Implementation 

The methods for calculating sensitivities and the example problems have been 
implemented using the general purpose finite element code, EAL [23]. EAL includes 
general language constructs for controlling execution flow as well as general and 
specific utilities for manipulating data stored as named entities in a database. It also 
allows procedures (called “runstreams”) to be defined and then explicitly executed. 
Most of the implementation was done using EAL runstreams. However, some parts 
of the implementation could not be conveniently done using runstreams and were 
coded as Fortran additions to EAL. The Fortran additions are described below. The 
runstreams for the algorithms and example problems are included and described 
below. 

Additions to EAL 

The transient response module in EAL version 312 solves the uncoupled form 
of equations (2.5) using the matrix series expansion method. A modification was 
made to allow equations (2.5) to be fully coupled. In the semi-analytical method, 
the right-hand side, pseudo loading of equations (4.6) can be easily formed using 
EAL. However, a slight modification to the transient response module was required 
to permit solution of equations (4.6) with this general form of loading. In addition, 
a special purpose module was added to EAL to perform the task of identifying the 
critical points on each response function. 

Runstream for Stepped Beam Example Problem 

The runstream for the stepped beam is included to illustrate how the sensi- 
tivity calculation runstreams are used. At the beginning of the runstream, the 
data set XFLG ADS indicates which subset of the possible design variables will be 
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considered in the sensitivity analysis. The data sets X ADS and XNAME ADS contain 
the initial values and the register names of all the design variables, respectively. 
Various parameters controlling the analysis and sensitivity calculations are defined 
in runstream data sets TR PARAMETERS, DXDV PARAMETERS, and BACK METHOD. The 
runstream data set MODEL defines the model in terms of the design variables in X 
ADS. It is called before the initial dynamic analysis and at least once for each design 
variable considered in the sensitivity analysis. The runstream data set DYNAM SOLN 
is called once to perform the dynamic analysis of the initial model. The runstream 
data set PLOT RESP illustrates the interface to a useful utility runstream TR PLOT 
for automatically generating plots of response quantities as a function of time. TR 
PLOT is called once for each class (e.g. accelerations) of response quantity to be 
plotted. The actual sensitivity analysis is performed by calling the runstream TR 
DXDV n where the n is associated with the particular sensitivity calculation method. 
*CM=120000 

t 

t RUNSTREAM FOR STEPED BEAM EXAMPLE 

* 

•XQT EITE 

! STST = SSP(4,S) $ GET SYSTEM TYPE 

•XQT U1 

*INF=7 

*CLIB=28 

•(ALL) ALL 

•XQT AUS 

TABLE ( NI=1 ,NJ=8 ,TYPE=0) : IFLG ADS 
J=1 : 1 
J=5 : 1 
J=6 : 1 
J=8 : 1 
•XQT U1 
•TI(X ADS) 

23. B 

22.0 

20.0 

18.0 
16.5 

* 

40.0 

80.0 

120.0 

160.0 
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*TI(XNANE ADS) 

HI : H2 : H3 : H4 : HS 
XL1 : XL2 : XL3 : XL4 
•(TR, PARAMETERS) 

QLIB-1 
MNAME-CEM 
NMODES-5 
DT=1.0E-S 
T2 * .3 
DRFORMAT=DIAG 
METHOD=MODES 
DXDV=0 
EIGEN=1 
PRINT® 1 
VLIB=1 
NCRIT=6 
C0NV=1 . E-10 
BLKSIZE=2000 
NTERMS=50 
* ( DXDV , PARAMETERS ) 

FDCH=1 . OE-B 
FDMCH=1.0E-6 
DXMD=FIXED 
*TI(BACK METHOD) 

2 $ DISP 
1 t VELOCITIES 
1 t ACCELERATIONS 

1 $ REACTIONS 

2 $ STRESSES 

•(MODEL) END 

! LEN=200 . 

NEPS = 3 
! NEL = NEPS*5 
! NNODE = NEL + 1 
•XQT ADS 

TABLE(NI=1,NJ=S) : XXI 

1=1 : J=1,B : 0. "XL1" "XL2" "XL3" "XL4" 
TABLE(NI=1 ,NJ=5) : XX2 

1=1 : J=1 ,5 : "XL1" "XL2" "XL3" "XL4" "LEN” 
D1 = SDM(XX2 -1.0 XXI) 

! RNEP = 1.0/NEPS 
DELX = UNION ("RNEP" Dl) 

•XqT TAB 

START "NNODE" 1346 
JLOC 

$1 0. 0. 0. 200. 0. 0. "NNODE" 

! J = 1 : ! X = 0.0 
! II = 1 i ! N1 = 5 
•LABEL 6 

! DELX = DS , 1 , "II" ,1(1 DELX ADS 1 1) 

! 12 = 1 : ! N2 = NEPS 
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•LABEL 8 
"J" "I" 0.0 0.0 
! J = J ♦ 1 : ! X = X ♦ DELI 
*JGZ,-1(N2,8) 

! II = II + 1 
*JGZ,-l(Ni ,6) 

"NNOD" "LEN” 0. 0. 

CON 1 

ZERO 1,2,6 : 1 
NATC 

1 30. +6 .3 .3 

BA 

RECT 1 1.20 "HI" 

RECT 2 1.10 "H2" 

RECT 3 1.00 "H3" 

RECT 4 .80 "H4" 

RECT 5 .85 "H6" 

RECT 6 1.0 20.0 
HREF 

1 1 2 1 1.0 
•XQT ELD 
E21 

! N = 5 
•1 = 1 
! N1 = 1 
•LABEL 20 
NSECT = "I" 

! N2 = N1 + 1 
"Nl" "N2" 1 "NEPS" 

! N1 = N1 ♦ NEPS 
! I = I + 1 
•JGZ,-1(N,20) 

• XQT E 
RESET G=386 . 

•XQT EKS 
•XQT TAN 
•XQT K 
•XQT H 
RESET G=386. 

•XQT ADS 
R = RIGID(l) 

DEFINE R6 = R ADS 1 1 6,6 
APPL FORC = PROD (-1.0 CEN R6) 
•END 

*(DTNAN,SOLN) END 
•DC ALL ( TR , VECTORS ) 

•XQT D1 
•TI(SEL DISP) 

"NNODE" 2 
•TI(SEL TELO) 

"NNODE" 2 



*TI(SEL 1CCE) 

"NNODE" 2 
*TI(SEL STRE) 

E21 1 1 SZ1 1 0 
*XQT AOS 

t DEFINE SOKE MODAL DAMPING 
TABLE(NI=1 ,NJ="NMODES") : DRAT 
1=1 : J=l,"NMODES" : .006 
•DC ALL ( SQUARE LOAD) RTIME=.18 RANG=10.0 
*DCALL(TR,MAIN) 

* END 

•(SQUARE, LOAD) END 
•xqT ADS 

! RT2 = RTIME/2 . 0 
! EPS = 0.0 
! RT2M = RT2 - EPS 
! RT2P = RT2 + EPS 
! RTM = RTIME - EPS 
! D2R = 3.1416926/180. 

! READ = RANG*D2R 
! AMP = 4 . 0*RRAD/RTIME/RTINE 
! MAMP = -AMP 
TABLE(NJ=6) : TIME 

J=1 ,6 : 0. "RT2M" "RT2P" "RTM” "RTIME" 10000.0 
TABLE(NJ=8) : CA 

J=1 ,6 : "AMP" "AMP" "MAMP" "RAMP" 0.0 0.0 
•END 

•(PLOT RESP) END 
•XQT DCD 

CHANGE 1 A ADS MASK MASK HIST CA 1 1 
•XQT ADS 
ALPHA : FTITLE 

1» HISTORY OF FORCE MULTIPLIER G(T) 

ALPHA : DTITLE 

1’ TIP DISPLACEMENT HISTORY FOR CANTILEVER BEAM 
ALPHA : TVTITLE 

1’ TIP VELOCITY HISTORY FOR CANTILEVER BEAM 
ALPHA : TATITLE 

1» TIP ACCELERATION HISTORY FOR CANTILEVER BEAM 
ALPHA : STITLE 

1> BENDING STRESS AT THE ROOT FOR THE CANTILEVER BEAM 
! TLIB = 16 
•XQT U1 

•(TRPLOT OPTIONS) 

! YNAME = ’CA 
! TITLE = ’FTITLE 
! ID = 1 
•DC ALL(TR , PLOT ) 

•XQT 01 

•(TRPLOT OPTIONS) 

! YNAME = ’DISP 



! IDJK = ’DISP 
! TITLE = ’DTITLE 
! ID = "NNODES" 

*DCALL(TR,PLOT) 

•XQT 01 

•(TRPLOT OPTIONS) 

! INANE * >TELO 
! IDJK = ’VELO 
! TITLE * ’TTTITLE 
! ID = "NNODES" 

•DC ALL (TR, PLOT) 

•xqT U1 

•(TRPLOT OPTIONS) 

! TNAME * ’ACCE 
! IDJK = ’ACCE 
! TITLE = ’TATITLE 
! ID = "NNODES" 

*DCALL(TJt,PLOT) 

•XQT 01 

•(TRPLOT OPTIONS) 

! INANE = >STRE 
! IDQ = ’STRE 
! TITLE = ’STITLE 
! ID = "NNODES" 

•DCALL(TR.PLOT) 

•END 

•RGI 

•DCALL (TR , PARAHETERS) 

•DC ALL ( SENS , DT0P ) 

•DCALL(MODEL) 

•DC ALL (DINAH, SOLN) 

$*DCALL ( PLOT , RESP ) 

•IFC'DXDI" NE 0): *DCALL(TR , DXDT , "DXDV" ) 
•ALL 

•IF("SIST" EQ CDC ): *PRT(ALL) 

•IFC'STST" EQ CNIX): *PRT(ALL) 
•DCALL(ALL) 

•XQT EXIT 


Runstreams for Sensitivity Methods 


Runstream TR MAIN 


This is the main runstream for performing the transient response analysis and 
is based on a similar runstream produced by EISI. It is used only for the transient 
analysis of the initial model and not in the sensitivity calculations. 

I 
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t (TR MAIN) - MAIN DRIYER FOR TRANSIENT RESPONSE ANALYSIS 

I 

* XQT 01 

* REGISTER ST0REC29 TR REGISTERS 11) 

* REGISTER RETR (29 TR REGISTERS 1 1) 

* RGI 

$ DEFAULT REGISTERS: 

qLIB =2 * SOURCE FOR EXCITATION, DESTINATION LIB FOR RESPONSE 

YLIB si $ SOURCE LIB FOR YIBR MODE AND YIBR EVAL DATASETS 

VSET * 1 $ USE "YLIB" YIBE MODE "VSET" "VCON" FOR THE RITZ 

YCON =1 $ FUNCTIONS 

KNAME = K $ STIFFNESS MATRIX 

HNAME = DEM $ MASS MATRIX 

DAMP = DAMP $ NAME OF SPAR FORMAT DAMPING MATRIX 

FSET =1 $ EXCITATION SET NUMBER 

NAME = AUS $ 2ND WORD OF TIME "NAME" AND CA "NAME" 

DT =0. $ SET TIKE INCREMENT 

T2 =0. $ FINAL INTEGRATION TIME 

DEFORM* DIAG t FORMAT FOR THE REDUCED MATRICES (DIAG, FULL, RITZ) 
DRMETH* 0 t TIME INTEGRATION METHOD (0=MSE) 

NTER * 50 $ SET NUMBER OF TERMS IN MATRIX SERIES EXPANSION 

NMODES= 0 I NUMBER OF MODES USED IN DYNAMIC ANALYSIS (DEFAULT=ALL) 

BLKSIZ= 0000 $ BLOCK SIZE FOR OUTPUT DATA SETS 

EIGEN =0 $ EIGENALUE ANALYSIS OF DAMPED SYSTEM 

PRINT = 0 $ PRINT FLAG FOR DTEX 

0PT=0 , PROC=MAIN, NERR=0 

* DC ALL, OPT (TR PARAMETERS) 

! ZERO = SSP(0,10) 

*IF("NMODES" EQ 0): ! NMODE=TOC ,NBLOCK("VLIB" YIBR MODE "YSET" "YCON") 

t 

t COMPUTE DATASETS REQUIRED FOR DR/DTEX, /TR1 , AND /BACK: 

$ 

* CALL (TR PREP) 

t 

t COMPUTE THE MODAL RESPONSE: 

I 

* XQT DRX 

* IF("DT" GT l.E-20) : *GOTO 20 

DTEX(INLI="QLIB",N2="NAME" , OUTL="QLIB" , EIGEN=»EIGEN" , > 

PRINT="PRINT" ) 

•GOTO 30 
•LABEL 20 

DTEX(INLI*"QLIB" ,N2="NAME" ,OUTL="QLIB" ,DTs"DT" ,> 
NTER="NTER",EIGEN="EIGEN",PRINT="PRINT") 

•LABEL 30 

TR1 ( INLIB="QLIB" , N2="NAME" , CASE=”FSET" , ALIB="QLIB"> 
QXLIB="QLIB" ,QX1LIB= "QLIB" ,QX2LIB= "QLIB" ,> 

T2="T2" ,LB="BLKSIZ") 

t 

t BACK TRANSFORMATION: 

t 
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♦XQT ADS 
! NBCK = 5 

TABLE(NI="NBCK», KJ=1, TTPE=4) : "QLIB" BACK LIST 
J=1 : DISP TELO ACCE &EAC STB.E 

t 

EXIT 

•1=1: ! N = NBCK 
* LABEL 50 

! BKMETH = DS,1,"I",1("QLIB» BACK METH 1 1) 

• NN = DS , "I" , 1 , 1 ( "QLIB” BACK LIST 1 1) 

! I EUR = TOC , IERR("QLIB" SEL "NN" MASK MASK) 
♦IF("IERR" EQ 0): ‘CALL (TR "NN" "BKMETH") 

! I = I ♦ 1 
*JGZ,-1(N,50) 

$ 

$ EXIT: 

• XQT U1 

• REGISTER RETRIEVE (29 TR REGISTERS 1 1) 

♦END 



Runstream TR PREP 


This runstream is used to define the reduced equations for the transient analysis 
and also prepare data for the back transformation phase. It is based on an EISI 
runstream of the same name. It is used only for the transient analysis of the initial 
model and not in the sensitivity calculations. 


(TR PREP) - PREPARATION OP REQUIRED DATASETS 


XQTC U1 
RGI 

PROC=PREP , NERR=0 , MOTI=MOTI, F0RC=F0RC 

CHECK FOR THE REQUIRED DATASETS AND DETERMINE THE TYPE OF EXCITATION: 
! TTPE=0 FOR APPLIED FORCE. !TTPE=1 FOR APPLIED DISPLACEMENT. 

!TTPE=2 

! IERR=TOC , IERR( "QLIB" APPL FORC "FSET" MASK) 

IFC’IERR" NE 0): *G0 TO 100 

!TTPE=0: !FNAME=’F0RC M APPLIED FORCE EXCITATION 
LABEL 100 

!IERR=TOC,IERR("QLIB" APPL MOTI "FSET" MASK) 

IFO'IERR" NE 0): *G0 TO 200 


IF ("TYPE" EQ 0): !NERR=1 M FORCE ft DISP SPECIFIED * NERR=1 

!TYPE=1: !FNAME=>MOTI M APPLIED DISPLACEMENT EXCITATION 
LABEL 200 

IF("TYPE" EQ 2) : !NERR=2 1$ NO EXCITATION SPECIFIED * NERR=2 

IFC'NERR" NE 0):*CALL (TR ERROR) 

!IERR=TOC,IERR("VLIB" VIBR MODE "VSET" "VCON"): !NERR=3 
IFC'IERR" NE 0) : *CALL (TR ERROR) I NERR=3 

! IERR=T0C , IERR( "VLIB" VIBR EVAL "VSET" "VCON"): !NERR=4 
IFC'IERR" NE 0) :*CALL (TR ERROR) # NERR=4 

! IERR=T0C , IERR( "QLIB" TIME "NAME" "FSET" MASK): INERR-5 
IFC'IERR" NE 0) :*CALL (TR ERROR) *NERR=6 

! IERR=T0C , IERR( "QLIB" CA "NAME" "FSET" MASK): !NERR=« 

IFC'IERR" NE 0) : ‘CALL (TR ERROR) $NERR=6 


: COMPUTE XTMX , XTKX, XTDX, ft ITF FOR DR/DTEX ft TR1 : 

! IERR=T0C , IERR( 1 INV "KNAME" "VCON" MASK) 
IF ("IERR" EQ 0) :*G0 TO 260 
XQT DRSI 

RESET K="KNAME" ,C0N="VC0N" 

LABEL 250 
XQT AUS 

OUTLIB="QLIB" : INLIB="QLIB" 
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DEFINE 

X="VLIB" 

VIBR 

MODE 

"TSET" 

DEFINE 

E=»TLIB" 

VIBR 

ETAL 

"VSET" 

DEFINE 

F="QLIB" 

APPL 

"FNAME" 

"FSET" 

DEFINE 

K= 1 

"KNAME 

H 


DEFINE 

H= 1 

"MNAME 

M 



* IF("TTPE" EQ 0): *G0 TO 300 

KS=PROD ( "KNAME" -1. F) : DEFINE F=KS 

• LABEL 300 

ITF "NAME" "FSET"s XTT(X,F) 

t 

! IDMD = TOC , IERR( "VLIB" DRAT MASK MASK MASK) 

*IF("IDMD" NE 0): *GOTO 400 

DEFINE D = "FLIB" DRAT 
OMEG = SQRT(E) 

DMPD = PRDD(2 . 0 D ONEG) 

•LABEL 400 

! IERR=TOC , IERR( "QLIB" XTMX "NAME" MASK MASK) 

* IF("IERR" EQ 0): *G0 TO 500 

$ 

•CALL(TR,REDM) 

• LABEL 500 

t 

$ COMPOTE ("QLIB" STAT DISP "FSET" "FCON"): 

t 

! IERR=TOC , IERR( "QLIB" STAT DISP "FSET" "VCON") 

• IF("IERR" EQ 0): *G0 TO 700 

* XQT SSOL 

RESET SET="FSET" , CON="VCON" , QLIB="QLIB" 

• LABEL 700 

* 

•LABEL 1000 
•END 
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Runstream TR DISP 


The two TR DISP runstreams do the back transformation for displacements. 

TR DISP 1 performs the back transformation using the mode displacement method 

and TR DISP 2 uses the mode acceleration method. This naming convention is 

used for the other runstreams that perform the back transformation operation for 

other response quantities. The TR DISP runstreams and companion runstreams for 

velocities, accelerations, and stresses perform the back transformation at all time 

steps. Accordingly they are used only in the dynamic analysis of the initial model 

and not in the sensitivity analysis. 

* 

t (TR DISP 1) - BACK TRANSFORMATION FOR DISPLACEMENTS 
$ MODE DISPLACEMENT METHOD 

$ 

•iqTC AOS 

OUTLIB="QLIB" : INLIB="C]LIB" 

DEFINE IDJK = "QLIB" SEL DISP 

DEFINE X = "VLIB" VIBR MODE "VSET" "VCON" 

TMAT YMOD=SVTRAN(IDJK , X) 

•XQT DRX 

BACK(LRZ="BLKSIZE") 

T = +1.0 "QLIB" TMAT VMOD : T = "QLIB" QX 

Z= "QLIB" HIST DISP 

EXT = "QLIB" EXT DISP "FSET" 

•END 

* 

$ (TR DISP 2) - BACK TRANSFORMATION FOR JOINT DISPLACEMENTS 

* USING MODE ACCELERATION METHOD 

* 

* 

$ TRANSIENT RESPONSE: BACK TRANSFORM FOR JOINT DISPLACEMENTS 
$ APPLICABLE FOR APPLIED FORCE OR DISPLACEMENT EXCITATION 

$ 

$ REGISTERS: QLIB, NAME, FSET, VCON 

* 

• XQT AUS 

! ID - TOC, IERR( "QLIB" XTDX MASK MASK MASK) 

0UTLIB="QLIB" : INLIB="QLIB" 

DEFINE E = "VLIB" VIBR EVAL 
ROMG = RECIP(E) 


DEFINE 

IDJK 

= "QLIB" 

SEL 

DISP 

DEFINE 

XS 

= "QLIB" 

STAT 

DISP "FSET" "VCON" 

DEFINE 

X 

= "VLIB" 

VIBR 

MODE 1 "VCON" 
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DEFINE DiCC = THAT DACC 
XQME = CBD(X.ROMG) 

THAT DACC = SVTRAN(IDJK.XOME) 

THAT DS ~ SVTRAN( IDJK , IS) 

*JNZ(ID , 30) 

$ DAMPING TEEM 

< NJ = TOC ,NJ("QLIB" XTDX MASK MASK MASK) 

! NINJ = TOC , NINJ ( "QLIB" XTDX MASK MASK MASK) 

•IFO'NJ" HE "NINJ"): *GOTO 10 
t MODAL DAMPING 

XOMD = CBD(XOME.XTDX) 

TMAT DVEL = SVTRAN ( IDJK , XOMD ) 

•GOTO 30 
•LABEL 10 
$ GENERAL DAMPING 

TMAT DVEL = RPROD(DACC, XTDX) 

•LABEL 30 
• XQT DRX 

BACK(LRZ="BLKSIZE") 

T = +1. "QLIB" TMAT DS : T = "QLIB" A "NAME" "FSET" 
•IFC'ID" EQ 0): T=-l . "QLIB" TMAT DVEL : T="QLIB" QI1 "NAME" "FSET" 
T = -1. "QLIB" TMAT DACC : T = "QLIB" QX2 "NAME" "FSET" 

Z = "QLIB" HIST DISP "FSET" 

EXT = "QLIB" EXT DISP "FSET" 

•END 


Runstream TR VELO 


* 

$ (TR VELO 1) - BACK TRANSFORMATION FOR VELOCITIES 
$ MODE DISPLACEMENT METHOD 

$ 

•XQTC AUS 

OUTLIB="QLIB": INLIB="QLIB" 

DEFINE IDJK = "QLIB" SEL VELO 

DEFINE X = "VLIB" VIBR MODE "VSET" "VCON" 
TMAT VVEL=SVTRAN(IDJK,X) 

•XQT DRX 

BACK (LRZ="BLKSIZE" ) 

T = +1.0 "QLIB" TMAT VVEL : T = "QLIB" QX1 

Z= "QLIB" HIST VELO 

EXT = "QLIB" EXT VELD "FSET" 


•END 



Runstream TR ACCE 


$ 

$ (TR ACCE 1) - BACK TRANSFORMATION FOR ACCELERATIONS 
t MODE DISPLACEMENT METHOD 

* 

*XQTC AUS 

OUTLIB="QLIB": INLIB="QLIB" 

DEFINE IDJK = "QLIB" SEL ACCE 

DEFINE X * "YLIB" YIBR MODE "VSET" "YCON" 
TMAT VACC=SVTRAN(IDJK , I) 

• XQT DRX 

BACK(LRZ="BLKSIZE") 

T = +1.0 "QLIB" TMAT YACC : T = "QLIB" QX2 
Z= "QLIB" HIST ACCE 
EXT s "QLIB" EXT 


•END 


ACCE "FSET" 



Runstream TR STRESS 


t 

I (TR STRESS 1) 

$ 

* 

t MODE DISPLACEMENT STRESS BACK TRANSFORMATION 

I 

•XQT ES 
RESET OPER-T 
IDQ= "QLIB" SEL STRESS 
0 = "VLIB" VIBR MODE 1 "VCON" l,”NMODE" 

T = ''QLIB" THAT VSTRE "FSET" 

•XQT DRX 

BACK(LRZ="BLKSIZE") 

T = +1. "QLIB" THAT VSTRE "FSET" : T = "QLIB" QX "NAME" "FSET" 
Z = "QLIB” HIST STRESS "FSET" 

EXT = "QLIB" EXT STRESS "FSET" 

• END 

* 

$ (TR STRESS 2) 

$ 

* 

$ TRANSIENT RESPONSE: BACK TRANSFORM FOR ELEMENT STRESSES 
t APPLICABLE FOR APPLIED FORCE OR DISPLACEMENT EXCITATION 

* 

t REGISTERS: QLIB, NAME, FSET, VCON 

$ 

* XQT AUS 

! ID = TOC, IERR( "QLIB" XTDX MASK MASK MASK) 

OUTLIB="QLIB" : INLIB="QLIB" 

DEFINE E = "VLIB" VIBR EVAL 
ROMG = RECIP(E) 

DEFINE IDJK = "QLIB" SEL 

DEFINE XS = "QLIB" STAT 

DEFINE X = "VLIB" VIBR 

XOME = CBD(X.ROHG) 

*JNZ(ID,30) 

* DAMPING TERM 

! NJ = TOC, NJ ("QLIB" XTDX MASK MASK MASK) 

! NINJ = TOC ,NINJ("QLIB" XTDX MASK MASK MASK) 

•IFC'NJ" NE "NINJ"): *GOTO 10 

* MODAL DAMPING 

XOMD = CBD( XOME, XTDX) 

•GOTO 30 
•LABEL 10 
I GENERAL DAMPING 

XOMD=CBR(XOME,XTDX) 


DISP 

DISP "FSET" "VCON" 
MODE 1 "VCON" 


•LABEL 30 



•XQT ES 


RESET OPER=T 

IDQ = "QLIB" SEL STRESS 

U= "QLIB" STAT DISP "FSET" "TCON" : T= "QLIB" THAT SF 
*IF("ID" EQ 0): U* "QLIB" XOMD : T= "QLIB" THAT SD 

H= "QLIB" IOKE : T= "QLIB" THAT SP 

* XQT DRX 

BACK (LRZ*"BLKSIZE" ) 

T » +1. "QLIB" THAT SF : Y * "QLIB" A "HAKE" "FSET" 
•IF("ID" EQ 0): T=-l. "QLIB" THAT SD : T * "QLIB" qXl "NAME" "FSET" 

T a -1. "QLIB" THAT SP : T = "QLIB" QI2 "NAME" "FSET" 

Z = "QLIB" HIST STRE "FSET" 

EXT * "QLIB" EXT STRE "FSET" 

• END 
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Runstream TR ERROR 


* 

$ (TR ERROR) 

* 

* 

* THIS PROCEDURE PRINTS FATAL ERROR MESSAGES FOR THE TR PROCS. 

t 

* IQT U3 

RP2: NUMBER OF F0RJUTS=10 

FORM 1 ’ (33H1*** TR FATAL ERROR: PROC, NERR= ,A4 P 1H,,I4) 

PRINT(l) "PROC" "NERR" 

* GO TO "PROC" 

$ 

* LABEL PREP 

* 

FORM 1 ’ (10X,47H .BOTH (APPL FORC FSET ) AND (APPL MOTI FSET )/ 

> 10X.46HARE SPECIFIED IN QLIB , ONLY ONE IS PERMITTED) 
FORM 2 ' ( 10X , 48HNEITHER (APPL FORC FSET ) NOR (APPL MOTI FSET )/ 

> 10X.20HIS PRESENT IN QLIB ) 

FORM 3>(10X,47HYIBR MODE YSET YCON IS NOT PRESENT IN YLIB ) 

FORM 4»(10X,47HVIBR EYAL VSET YCON IS NOT PRESENT IN YLIB ) 

FORM S’(10X,43HTIME NAME FSET IS NOT PRESENT IN QLIB ) 

FORM «>(10X,43HCA NAME FSET IS NOT PRESENT IN QLIB ) 

PRINT ("NERR") 

$ 

* GO TO FINIS 

$ 

* LABEL FINIS 

* XQT U1 

* SHOW 

* XQT DCU 

TOC 1: TOC "QLIB": TOC "YLIB" 

* XQT EXIT 
•END 
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Runstream TR RITZ 


This runstream calculates RWL vectors following equations (2.25) through 
(2.29). Then the reduced is system is optionally uncoupled by solving the reduced 
order eigenproblem. This runstream is substantially based on one written at EISI. 

$ 

$ (TR RITZ) 

* 

* xqT ui 

* 

$** REGISTERS: HAXHZ, AFLIB, VLIB, KNAME, NMD, SCALE 

$ 

* REGISTER STORE (29 REGISTER HORSE 1 1) 

« REGISTER RETR (29 REGISTER HOUSE 1 1) 

!IERR=TOC IERR(1 INY K MASK MASK) 

*ONLINE=0 

* IFO'IERR" EQ 0): *G0 TO 109 

* xqT DRSI 

* LABEL 109 

* xqT SSOL 

RESET qLIB="AFLIB" 

* xqT AUS 

0UTLIB=10: INLIB=10 
DEFINE M=1 "MNAME" 

I 

$ SCALE THE FIRST VECTOR 

t 

DEFINE X="AFLIB" STAT DISP 1 MASK 1 
MX=PR0D(M,I) 

XTMX=XTT(X,MX) 

RECI=RECI(XTMX) 

SCAL=SqRT (REC I ) 

11 RITZ YECT=CBD(X,SCAL) 

DEFI RITZ=11 RITZ VECT 

12 MX=PROD(M,RITZ) 

!NSET=T0C NBL0CKS(" AFLIB" STAT DISP 1 MASK) 

!N=NSET-1 : !N1=1 : !N2=0 

* IF ("N" Eq 0): *G0 TO 104 

$ H-ORTHONORMALIZE VECTORS 2 THROUGH NSET 

* LABEL 10S 

!N1=N1+1: !N2=N2+1 

DEFI U="AFLIB" STAT DISP 1 MASK "Nl" 

DEFI MX-12 MX MASK MASK MASK 1 "N2" 

XTMU=XTT(MX,U) 

DEFI X=ll RITZ VECT MASK MASK 1 "N2" 

A=CBR(X , XTMU) 
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U1=SUM(U,-1. A) 

MU=PR0D(M,U1) 

utmu=xttd(ui,mu) 

reci=reci(utmu) 

SCAL=SQRT(EECI) 

VECT=CBD (U1 , SCAL) 

11 EITZ VECT*UNION,U(VECT) 

TEMP = PROD ( M , VECT ) 

12 MX =UNION,U(TEMP) 

* JSZ -1 (N 106) 

* LABEL 104 

* XQT DCU 

ERASE 10 

! NDO=MMD-NSET: !SET=NSET*1: !SET1=HSET 

* IF ("NDO" LE 0): *G0 TO 1002 

* XQT AUS 

TABLE(NJ=1) : 13 SCALE 
J=l: "SCALE" 

$ 

t GENERATE REMAINING VECTORS ORTHONORMAL TO STATIC SOLUTION RITZ VECTORS 

I 

* LABEL 1000 

* XQT AUS 

0UTLIB=10: INLIB=10 
DEFINE M=1 "MNAME" 

DEFINE X=ll RITZ VECT MASK MASK "SET1" 

TEMP=PROD(M,X) 

NORM= NORM ( TEMP ) 

DEFI SCALE= 13 SCALE 
APPL FORC=CBD (NORM, SCALE) 

* IF(”SET1" EQ "NSET"): *G0 TO 106 
12 MX=UNION,U(TEMP) 

* LABEL 106 

* XQT SSOL 

RESET QLIBslO 

* XQT AUS 

OUTLIBslO: INLIB=10 
DEFI U=STAT DISP 

DEFI MX=12 MX MASK MASK MASK 1 "SET1" 

XTMU=XTT ( MX , U ) 

DEFI X=ll RITZ VECT MASK MASK 1 "SET1" 

A=CBR(X,XTMU) 

U1=SUM(U,-1. A) 

DEFI M=1 "MNAME" 

MU=PR0D(M,U1) 

UTMU=XTTD(U1 ,MU) 

RECI=RECI(UTMU) 

SCAL=SQRT(RECI) 

VECT=CBD(U1 , SCAL) 

11 RITZ VECT =UNION , U ( VECT ) 

!SET1=SET: !SET=SET+1 



• IQT DCtr 

ERASE 10 

• JGZ, -KUDO, 1000) 

• LABEL 1002 

* IF ( "DRFORMAT" EQ DIAG): *GOTO 10020 

• IQT ADS 

DEFI 1=11 RITZ VECT 

"VLIB" VIBR MODE 1 1*DNI0N(X) 

TABLE(NI=1 ,NJ*"NMD") s VIBR EVAL 
•RETDRN 
•LABEL 10020 

• IQT ADS 

0DTLIB=10: INLIB=10 

DEFINE K=1 K: DEFI M=1 "MNAME" 

DEFINE X=ll RITZ VECT 

IJC0DE=10000 

!NMODE=NMD 

KX=PROD(K,X): STN K 10000 "NMODE" = XTTS(X.KX) 
MX=PROD(M,X) : STN H 10000 "NMODE" = ITTS(X,MX) 
! ZER0=NM0DE-1 

• JZ (ZERO, 1003) 

•xgT DCD 

TOC 10 

• XQT STRP 

RESET S0DRCE=10, DEST=10 , FRQ2="MAXHZ" 

• JGZ (ZERO, 1004) 

• LABEL 1003 

• XQT ADS 

0DTLIB=10: INLIB=10 

!K=DS 2 1 1(10 STN K MASK MASK) 

!M=DS 2 1 1(10 STN M MASK MASK) 

!EVAL=K/M 

TABLE(NI=1 ,NJ=1) : STS EVEC: J=l: 1.0 
TABLE(NI=1,NJ=1): STS EVAL: J=l: "EVAL" 

• LABEL 1004 

• XQT ADS 

0DTLIB=10: INLIB=10 
DEFINE E=STS EVEC 
DEFI X=ll RITZ VECT 
X ORTH 1 1=CBR(X,E) 

DEFINE I=X ORTH 1 1 

"VLIB" VIBR MODE 1 1=DNI0N(X) 

DEFINE E=STS EVAL 

"VLIB" VIBR EVAL 1 1=DNI0N(E) 

"VLIB" VIBR HZ 1 1=SQRT(. 0253303 E) 

* ONLINE* 1 

• XQT DCD 

PRINT "VLIB" VIBR HZ 1 1 

• XQT N1 

• REGISTER RETR (29 REGISTER HODSE 1 1) 

•END 


152 



Runstream TR REDM 


This is a utility runstream for generating the reduced equations given a set 
of basis vectors. Depending on the input register DRFORMAT the equations can be 
coupled or uncoupled. If the equations are uncoupled, it is assumed that 
is the identity matrix and f* T K# is a diagonal matrix with the eigenvalues along 
the diagonal. 

* 

* (TR REDM) - FORM REDUCED K AND N MATRICES FOR TRANSIENT RESP. 

* 

t 

I REGISTERS: 

$ DRFO = 'FULL, >RITZ, OR 'DIAG 

$ NMODE = NUMBER OF MODES 

$ MNAME = MASS MATRIX NAME 

I VLIB = LIBRARY FOR VIBRATIONAL MODES AND FREQS 

t QLIB = DESTINATION LIBRARY FOR MATRICES 

t 

•XQTC AUS 
OUTLIB»"QLIB" 

DEFINE X * "VLIB" VIBR MODE 1 1 1,"NM0DES" 

DEFINE E * "VLIB" VIBR EVAL 1 1 

DEFINE DAMP = 1 DAMP SPAR 
DEFINE DMPD = 1 DMPD 

! IDSP = TOC ,IERR(1 DAMP SPAR MASK MASK) 

! IDMD = TOC, IERR( 1 DMPD MASK MASK MASK) 

! DRFO 

• IFO'DRFO" NE FULL): *G0T0 100 

t 

$ FULL MATRICES, X IS A SET OF EIGENVECTORS 

$ 

! N = NMODES 
! I = 1 

TABLE(NI="NMODES" ,NJ="NMODES") : XTMX 
•LABEL 10 

I=»I» ; J=»I» : 1.0 
! I = I + 1 
•JGZ ,-l(N , 10) 

! N = NMODES 
! I = 1 

TABLE (NI= "NMODES" , N J = "NMODES " ) : XTKX 
•LABEL 20 

! K = DS, "I", 1,1 ("VLIB" VIBR EVAL 1 1) 

I="I" : J="I" : "K" 

II = 1 + 1 
*JGZ,-1(N,20) 
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•IFC'IDSP" HE 0): *G0T0 8B 
0UTL=22 : INLI=22 
DX x PR0D(DAMP,X) 

"QLIB" XTDX = XTT(X.DX) 

* LABEL 86 

*IF("IDMD" HE 0): *RETURN 
! H x NHODES : ! I a 1 

• IFC'IDSP" HE 0): TABLE (HI="HMODES" , HJ = "HMODES" ) : "QLIB" XTDI 
•IFC'IDSP" Eq 0): TABLE, U : "QLIB" XTDX 
•LABEL 00 

! D x DS,"I",l,l("VLIB" DHPD MASK MASK MASK) 

. j_«jn . iiqii 

! I x I 4 1 

•JGZ,-1(H,90) 

•RETURH 
•LABEL 100 

t 

$ FULL REDUCED MATRICES (X HOT EIGENVECTORS) 

* 

•IFC'DRFO" NE RITZ): *GOTO 200 
0UTLIB=22 : INLIBx22 
DEFINE K x 1 K SPAR 
DEFINE M x 1 "MNAME" 

KX = PROD(K,X) 

MX = PROD(M.X) 

"QLIB" XTKX x XTT(X,KX) 

"QLIB" ITMX x ITT(X.MX) 

•IFC'IDSP" NE 0): *GOTO 130 
0UTL=22 : INLI=22 
DX x PROD (DAMP, I) 

"QLIB" XTDX = XTT(X,DX) 

•LABEL 130 

•IF("IDMD" NE 0): •RETURN 
! N x NMODES : ! I = 1 

•IFC'IDSP" NE 0): TABLE (NI="NMODES" ,HJ="HMODES") : "QLIB" XTDX 
•IFC’IDSP" EQ 0): TABLE, U : "QLIB" XTDX 
•LABEL 140 

! D x DS,"I",1,1("VLIB" DMPD MASK MASK MASK) 

I="I" : Jx"I" : "D" 

! I = I + 1 
*JGZ,-1(N,140) 

•RETURN 
•LABEL 200 
8 

$ SIMPLE DIAGONAL CASE (X EIGENVECTORS) 

$ 

TABLE(NIxl,NJx"NMODES") : XTKX : TRAN(SOUR=E) 

TABLE(NI=1 ,NJ="NMODES") : XTMX : Jxl, "NMODES" : 1.0 
•IFC'IDSP" NE 0): *GOTO 210 
0UTL=22 : INLIx22 
DX = PROD (DAMP ,X) 



"QLIB" ITDX = XTTD(X,DX) 

* LABEL 210 

*IF("IDMD" HE 0): *RETURN 
! M = NHODES : ! I * 1 

*IF("IDSP" NE 0): TABLE(HI=1 ,NJ="NMODES") : ''QLIB" ITDX 
*IF("IDSP" Eq 0): TABLE, 0 : "QLIB" XTDX 
•LABEL 220 

! D = DS , "I " , 1 , 1 ( " VLIB" DHPD MASK MASK MASK) 

1=1 : 

! I = I ♦ 1 
*JGZ,-l(N,220) 

•RETURN 

•END 



Runstream TR PLOT 


This is a utility runstream for producing plots of response quantities as a func- 
tion of time. Its use is demonstrated in the stepped beam example runstream. 

t 

t (TR PLOT) 

t 

♦XQT 01 

# 

$ PLOTS TRANSIENT , TIME BISTORT DATA PRODUCED BT DR/TRI 

$ 

♦REGISTER EXCEPTIONS TLIB 
♦REGISTER STORE ("TLIB" TR REG 1 1) 

♦REGISTER RETREIVEC'TLIB" TR REG 1 1) 

* DEFAULT REGISTER ASSIGNMENTS 
! INLIB=1 

* IDJK = ’NONE 
! IDQ = ’NONE 

! TNAME = ’DISP 
! N3 = 1 

! TITLE = ’TITLE 
! ID = 1 
! 0PT=0 

♦DATA , OPTfTRPLOT OPTIONS) 

I 

! NS1 = TOC , NI ( " INLIB" HIST "TNAME"" N3” MASK) 

! NW1 = TOC ,NWDS("INLIB" HIST "TNAME" "N3" MASK) 

! NBLK = TOC ,NBLOCKS ("INLIB" HIST "TNAME" "N3" MASK) 

! NJBL = TOC, NJ(" INLIB" HIST "TNAME" "N3" MASK) 

! NS1 : .' NW1 : ! NBLK : ! NJBL 
! NSTE=NV1/NS1 
! NPPT=NSTE 

! DT=DS ,1,1 , 1 ("INLIB" DT MASK MASK MASK) $ TIME STEP 
! TIDQ = TOC ,IERR("INLIB" SEL "IDQ" MASK MASK) 

! TIDJ = TOC, IERR( "INLIB" SEL "IDJK" MASK MASK) 

! TTIT = TOC, IERR( "INLIB" "TITLE" MASK MASK MASK) 

*0NLINE=0 
♦xqT AUS 

TABLE(NI=1,NJ="NPPT") : "TLIB" XTAB 
DDATA="DT" 

J=1,"NPPT" : 0.0 

t 

t LOOP OYER ALL RESPONSE QUANTITIES 

* 

! NJLS = 1-NBLK*NJBL ♦ NSTE t NJ OF LAST BLOCK 
! KBLK = NBLK - 1 
! DBLS = KBLK*NJBL 
! JBLK = KBLK 
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! NSM1 = NS1 - 1 
! SBASE = 0 
! NJLS 

! 11=1 : ! N1=NS1 

* LABEL LI 

DEFINE T = "INLIB" HIST "YNAHE" "N3" 1 1,"JBLK" 

TABLE (NJ="NSTE") : "TLIB" TTAB ADS "II" 

*JZ(KBLK ,L2) 

TRAN(S0UR=T,SBAS="SBAS",SSKIP="NSH1",ILIM=1,JLIM="NJBL") 

•LABEL L2 

DEFINE T = "INLIB" HIST "TNAME" "N3" 1 "NBLK" , "NBLK" 

TABLE, D : "TLIB" TTAB AUS "II" 

TEAN(S0D&=T,SBAS="SBAS",SSKIP="NSN1",ILIM=1,JLIH="NJLS",DBAS="DBLS") 
! SBASE = SBASE ♦ 1 
! 11 = 11+1 
*JGZ,-1(N1,L1) 

$ 

t GENERATE AN X,Y PLOT FOR EACH RESPONSE QUANTITY 

* 

•ONLINEs 1 
•iqT PXY 

RESET DEVICE=META 
RESET NDEV=4014 
FONT XNDN=1 : FONT YNUH=1 
FONT ILAB=1 : FONT TLAB=1 
FONT TEIT=1 
X = "TLIB" XTAB 
XLABEL’ TINE (SECONDS) 

• JNZ (TID J ,110) 

TLFORMAT , 72> 

’ (4H J =,I2,9H JOINT = ,I6,8H COHP = ,I2,8H HIST = ,A4,6H ID = ,18) 
•LABEL 110 
•JNZ(TIDQ , 120) 

TLFORMAT, 72> 

>(4H J= , 12 , IX, A4,8H GRP= ,I2,6H IND= ,I6,7H COMP= ,A4,5H ID= ,16) 
•LABEL 120 
XAXIS=3,5 , 10 
YAXIS=4 , S ,10 
TP0S=0,0 

! 11=1 : ! N1=NS1 
•JNZ(TTIT,L3) 

TEXT = "TITLE" 

•LABEL L3 

ADVANCE 

BOUNDARIES = .01 .99 .04 .1 
*JNZ(TTIT,L4) 

PLOT TEXT 

•LABEL L4 

BOUND ARIES= .01 .99 .IS .85 
Y = "TLIB" TTAB AUS "II" 

•JNZ(TIDJ , 210) 
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! JOINT = DS , 1 , "I 1" , 1 ("INLIB" SEL "IDJK" MASK MASK) 

! COMP = DS,2,"I1",1("INLIB" SEL "IDJK" MASK MASK) 

TLABEL "II" "JOINT" "COMP" "TNAME" "ID" 

•LABEL 210 
•JNZ(TIDQ,220) 

! ENAME = DS,1, "II", 1 ("INLIB" SEL "IDQ" MASK MASK) 

! EGRP = DS,2,"I1",1("INLIB" SEL "IDQ" MASK MASK) 

! EINDX = DS,3,"I1",1("INLIB" SEL "IDQ" MASK MASK) 

! ECOMP = DS , 4 , "I 1" , 1 ( "INLIB" SEL "IDQ" MASK MASK) 

TLABEL "II" "ENAME" "EGRP" "EINDX" "ECOMP" "ID" 
•LABEL 220 
INIT 

PLOT CURT 
! 11 * 11+1 
•JGZ ,-l(Nl ,L3) 

•IQT 01 

•REGISTER RETRIEVE ("TLIB" TR REG 1 1) 

•FREE "TLIB" 

•RETORN 

•END 



Runstream TR VECTORS 


This runstream generates basis vectors by calling the system eigensolver, or 
calling runstream TR RITZ, or using the static mode method, or by other experi- 
mental techniques. 

$ 

* (TR, VECTORS) - COMPUTE VECTORS FOR USE IN DTNAMIC ANALYSIS 

* 

*XQT AUS 

R x RIGID(l) 

CR = PROD ( "MNAME" , R) 

Z = NDDF.l(CR) 

! NDDF = DS, 1,1, 1(1 Z AUS 1 1) 

! NDDF 

•IFC’HETHOD" NE MODE): *G0T0 100 
«XqT E4 

RESET NM0DES="NM0DES" 

RESET M="MNAME" 

RESET NDDF="NDDF" 

RESET C0NV="C0NV" 

•DC ALL, OPT (E4 PARAMETERS) 

•GOTO 200 

* 

•LABEL 100 

• IFC’HETHOD" NE RITZ): *G0T0 106 
! MAXHZ = 1 . 03E+10 
! NMD = NHODES 
! SCALE =1.0 
! AFLIB = 1 
! VLIB = 1 
*DCALL(TR RITZ ) 

•GOTO 200 

$ 

•LABEL 106 

•IFC’METHOD" NE OLD): *G0T0 110 
•XqT DCU 

COPT 3 1 VIBR MODE 
COPT 3 1 VIBR EVAL 
•GOTO 200 
•LABEL 110 

•IFC'HETHOD" NE ONES): *GOTO 116 
*DCALL(TEST NEB3) 

•GOTO 200 
•LABEL 116 

•IFC’METHOD" NE STAT): *GOTO 120 
•XqT E4 

RESET NMODES="NMODES" 

RESET H="MNAHE" 
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RESET NDDF="NDDF" 

RESET CONV="CONV" 

•DC ALL, OPT (E4 PARAMETERS) 

•XQT DRSI 
•XQT SSOL 

$ ORTHOGOGONALIZE STATIC SOLDTION AND APPEND TO SET OF MODE SHAPES 
•DC ALL (TR, GRAM) 

$ HAKE THE VECTORS ORTHOGONAL WITH RESPECT TO BOTH X AND M 
*DCALL(TR DIAG) 

•xqT VPRT 
PRINT SI ADS 
•XQT DCD 
TOC 1 

PRINT 1 VIBR EVAL 
•LABEL 120 

•IFC'METHOD" NE DMOT) : *G0 TO 200 
•xqT ADS 
Z = NDDP, 1 , 2(CR) 

•XqT DRSI 
RESET C0N=2 
•XqT ADS 
DDF = 1 
SSPREP(K,2) 

•XqT DCD 

CHANGE 1 BNF MASK MASK MASK VIBR MODE 1 1 
$ DNCODPLE THE SYSTEM 
•DCALL(TR.DIAG) 
t END OF METHODS OPTIONS 
•LABEL 200 
• END 



Runstream TR GRAM 


This is a utility runstream for performing a Gram-Schmidt orthogonalization 
of a set of vectors. 


(TR.GRAM) - PERFORM GRAM-SCHMIDT PROCESS TO M-ORTHOGONALIZE 
STAT DISP WITH RESPECT TO YIBR MODE AND THEN 
REPLACE THE LAST YIBR MODE WITH THE NEW VECTOR 


*XQT ADS 

! NMM1 = NMODES - 1 
DEFINE X = YIBR MODE 1 1 1,"NMM1" 

DEFINE S = STAT DISP 1 1 
DEFINE M = "MNAME" 

INLIB=10 : 0UTLIB=10 

$ ORTHOGONALIZE THE STATIC SOLUTION WITH RESPECT TO THE MODE SHAPES 
MX = PROD (M, I) 

XTMS = XTT(MX,S) 

A = CBR(X.XTMS) 

SI = SUM(S, -1.0 A) 

$ NOW SCALE THE VECTOR 
MS = PR0D(M,S1) 

STMS s ITY(Sl.MS) 

• STMS = DS, 1,1, 1(10 STMS AUS 11):! STMS * STMS**.S : ! STMS = 1.0/STMS 
TEMP MODE 11= UNIONO'STMS" S1,X) 

*XQT DCD 

CHANGE 10 TEMP MODE 1 1 VIBR MODE 1 1 
COPT 10 1 YIBR MODE 1 1 
•DELETE 10 
•END 
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Runstream SENS DVTJP 


This is a utility runstream for updating the design variable registers based on 
the data sets X ADS and XNAME ADS. It is always called immediately before calling 
model so the current values of the design variables are available for use. 

t 

* (SENS DVUP) - UPDATE DESIGN VAEIABLE REGISTERS FROM DATASET 

t 

•xqT ui 

! N = T0C,NJ(1 INANE ADS 1 1) 

! I = 1 
•LABEL 10 

! RNAHE = DS , 1 ,"I" ,1(1 INANE ADS 1 1) 

! RVAL = DS , 1 , "I" ,1(1 I ADS 1 1) 

! "RNANE" = "RVAL" 

! "RNANE" 

! I = I ♦ 1 
•JGZ ,-l(N, 10) 

•END 
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Runstream TR DXDV 1 


The TR DXDV n runstreams implement the different sensitivity methods. The 
structure of all these runstreams is similar. In each case there is a loop over the 
designated design variables and sensitivities of the required response quantities at 
the set of critical points are calculated. Within this loop there is at least one call 
to runstream MODEL to form a perturbed design, a call to form a set of new reduced 
equations, and a call to processor DRX to integrate the reduced equations in time. 
Runstream TR DXDV 1 implements the forward difference method with either fixed 
or updated basis vectors. 

* 

* (TR DXDV 1) - CALCULATES DERIVATIVES OF TRANSIENT RESPONSE 

$ USING THE FORVARD DIFFERENCE OPERATOR AND 

$ EITHER FIXED OR UPDATED MODES 

t UPDATE HISTORY 

t 6/28/88 WHG - MODIFIED FOR VELO, ACCE, STRESSES 

I- 

•XQT U1 

* REGISTER STORE (1 DXDV REGISTERS 1 1) 

•REGISTER RETRIEVE(1 DXDV REGISTERS 1 1) 

•RGI 

FDCH = .001 
FDMCH = .0001 
XLIB = 5 
OPT = 0 
RLIB = 14 
DRHETH0D=0 
DXMD=UPDATED 

•DC ALL, OPT (DXDV PARAMETERS) 

* SHOU 

•DCALL(TR.DPREP) 

$ 

t LOOP OVER ALL DESIGN VARIABLES 

* 

! NDV = TOC , NJ ( 1 X ADS 1 1) 

! NCNT = NDV 
! IDV = 1 
•LABEL 10 
•XQT U1 

! IFLG = DS , 1 ,"IDV" ,1(1 XFLG ADS 1 1) 

•JZ(IFLG,100) 

•LIBS "XLIB" 2 3 4 1 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
•XQT DCU 

COPT "XLIB" 1 XNAME ADS 1 1 
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•XQT U3 
RP2 

FORMAT 1>(1H1,20X,41HBEGINNING SEMSITITITT CALCULATION. DT NO., 13) 
PRINT(l) "IDT" 

*XQT AUS 

DEFINE X « "XLIB" X ADS 1 1 
TABLE(NI=1,NJ="NDT") : X ADS 1 1 
J=1 , "NDT" : 1.0 

J*"IDT FDCH" 

TRAN ( SOURCE 2 ! , OPERATION=NULT ) 

$ CHECK FOR TOO SMALL A STEP 
! X = DS , 1, "IDT", 1 ("XLIB" X ADS 1 1) 

! DX » FDCH*I 

*IF("DX" GT "FDCHM") : *GOTO 20 
! DX = FDMCH 
! X = X + FDMCH 
TABLE, U : I ADS 1 1 
OPER = XSUM 
J="IDT" : "X" 

•LABEL 20 
•CALL(SENS.DTUP) 

$ FORM PERTURBED MODEL 

*0NLINE=0 

•DCALL(MODEL) 

*0NLINE=1 
•XQT DCU 

COPT "XLIB" 1 TIME 
COPT "XLIB" 1 CA 
COPT "XLIB" 1 DMPD 
•IFC'DXMD" NE FIXE): *GOTO 30 
COPT "XLIB" 1 TIBR MODE 
•DCALL(TR.DIAG) 

•GO TO 40 
•LABEL 30 

•IFC’DXMD" NE UPDA): *GOTO 40 
•DC ALL ( TR , TECTORS ) 

•LABEL 40 
•XqT AUS 

DEFINE X = TIBR MODE 1 1 1,"NM0DES" 

DEFINE F = APPL FORC 1 
XTF = XTT(I,F) 

•DCALL(TR.REDM) TLIB=1 
•XQT DRI 

DTEI(DT="DT",METHOD="DRMETHOD",NTERMS="NTERMS") 

TR1 (QXLIB=1 , QX1L=1 , QX2L=1 ,T2="T2" ,LB="BLKSIZE") 

•DCALL(TR,DBACK 1) 

$ 

t COMPUTE DERITATITES USING FORWARD DIFFERENCE OPERATOR 

$ 

! OTDI = 1.0/DX 
! MOTD = - OTDX 
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! I = 1 : ! N = HBCK 

•LIBS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
•LABEL 80 
•XqT AUS 

! NM = DS, "I", 1,1(1 BACK LIST 1 1) 

! I ERR = T0C,IERR(1 SEL "NM" MASK MASK) 

•JNZ(IEBB,90) 

DEFINE CPI = "XLIB" CRPT "NM" 

DEFINE CPO = CRPT "NM” 

DXDT "NM" "IDT" = SUM("OTDX" CPI "MOVD" CPO) 

•XqT DC0 

PRINT 1 DXDT "NM" "IDT" 

•LABEL 90 
! I * I ♦ 1 
•JGZ,-1(N,80) 

ERASE "XLIB" 

•LABEL 100 
! IDT = IDT ♦ 1 
*JGZ,-1(NCNT,10) 

* 

•xqT Di 

•REGISTER RETRIETE(1 DXDT REGISTERS 1 1) 

•RETURN 

•END 
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Runstream TR DXDV 3 


This runstream implements the fixed mode semi-analytical sensitivity method. 
Depending on the call to runstream TR DBACK either the mode displacement or 
mode acceleration method is used to recover the physical sensitivities. 

$ 

* (TR DXDV 3) - CALCULATES DERIVATIVES OF TRANSIENT RESPONSE 
$ SEMI -ANALYTICALLY 

t 

*XQT 01 

* REGISTER ST0RE(1 DXDV REGISTERS 1 1) 

♦REGISTER RETRIEVE ( 1 DXDV REGISTERS 1 1) 

•RGI 

FDCH = .001 
FDMCH = .0001 
XLIB = 5 
OPT = 0 
RLIB = 14 
DRMETH0D=0 

•DC ALL, OPT (DXDV PARAMETERS) 

•SHOW 

$ 

I INITIALIZATION FOR DERIVATIVE CALCULATIONS 

$ 

•DCALL(TR.DPREP) 

•IQT U1 

t 

$ LOOP OVER ALL DESIGN VARIABLES 

$ 

! NDV = T0C,NJ(1 X ADS 1 1) 

! NCNT = NDV 
! IDV = 1 
•LABEL 10 

•LIBS "XLIB" 2 3 4 1 6 7 8 9 10 11 12 13 14 IS 16 17 18 19 20 
•XQT U1 

! IFLG = DS, 1, "IDV", 1( "XLIB" XFLG ADS 1 1) 

*JZ(IFLG , 100) 

•XqT DCU 

COPT "XLIB" 1 XNAME ADS 1 1 
•XQT U3 
RP2 

FORMAT 1 ’ (1H1 , 20X ,41HBEGINNING SENSITIVITY CALCULATION. DV NO. ,13) 

PRINT(l) "IDV" 

•XQT AUS 

DEFINE I = "XLIB" X ADS 1 1 
TABLE(NI=1,NJ="NDV") : X ADS 1 1 
J=1,"NDV" : 1.0 
J="IDV" : "FDCH" 


166 



TRAN(SOURCE=X, OPERATION=MULT ) 

I CHECK FOR TOO SHALL A STEP 
! I = DS , 1 , "IDT" , 1 ( "XLIB" X ADS 1 1) 

! DX = FDCH*X 

*IF("DX" GT "FDCHM") : *GOTO 20 
! DX - FDMCH 
! X * X ♦ FDMCH 
TABLE ( H : X ADS 1 1 
OPER = XS0H 
J="IDV" : "X" 

•LABEL 20 
•CALL(SENS,DVOP) 

$ FORM PERTURBED MODEL 

•0NLINE=0 

•DC ALL ( MODEL) 

*0HLINE=1 
•XQT ADS 

DEFINE I = "XLIB" VIBR MODE 1 1 1,"NH0DES" 

DEFINE KO = "XLIB" K SPAR 

DEFINE K1 = K SPAR 

DEFINE MO = "XLIB" "MNAME" 

DEFINE Ml a "MNAME" 

DEFINE DO = "XLIB" DAMP SPAR 

DEFINE D1 = DAMP SPAR 

DEFINE FO a "XLIB" APPL FORC 1 

DEFINE FI a APPL FORC 1 

! IDSP a TOC , IERR(1 DAMP SPAR MASK MASK) 

! OVDX a 1.0/DX 
! MOVD a -OVDI 

DKDV = SOM("OVDX" K1 "MOVD" KO) 

DMDV a SUM("OVDI" Ml "MOVD" MO) 

DFDV = SOH("OVDX" FI "MOVD" FO) 

•IFC’IDSP" EQ 0): DDDV a SUH("OVDX" D1 "MOVD" DO) 

DKX a PROD (DKDV, X) 

DMX a PROD (DMDV, X) 

•IF("IDSP" EQ 0): DDX = PROD(DDDV.X) 

RDKI = ITT (X, DKX) 

RDMI a ITT(X,DMX) 

•IFC'IDSP" EQ 0): RDDX a ITT(I.DDX) 

XTF AOS a ITT(X,DFDV) 

•XQT DCO 

COPT "XLIB" 1 TIME 
COPT "XLIB" 1 CA 
COPT "XLIB" 1 DT AUS 
COPT "XLIB" 1 DTEX AOS 

COPT "XLIB" 1 DCON AOS t CONSTANTS FOR NEWMARK METHOD 

I 

$ FORM THE RIGHT-HAND-SIDE PSEODO LOAD VECTOR 

* 

•XQT DRX 

BACK(LRZ="BLKSIZE" ,PRINT=0) 
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T = -1.0 MUX : T = "XLIB” QX2 XUS 
•IFO'IDSP" EQ 0): T = -1.0 MDX : T = "XLIB" QX1 AOS 
T = -1.0 MKX : T = "XLIB" QX AOS 
Z : FH AOS 
•XQT DEX 

TRl(QXLIB=l,qilLIB=l,QI2LIB=l,T2="T2",FHLIB=l,LB="BLKSIZE") 

$ 

$ BACK TRANSFORM FOR NECESSARY SENSITIVITIES OF PHYSICAL 
t QOANTITIES 

* 

•DC ALL (TR , DB ACK ,4) 

•XqTC DCO 
ERASE 1 

* LAB EL 100 

! IDV = IDV ♦ 1 
*JGZ,-1(NCNT,10) 

* 

•LIBS 1 2 3 4 S 6 7 6 9 10 11 12 13 14 IS 16 17 18 19 20 

•xqT 01 

•REGISTER RETRIEVEd OXDV REGISTERS 11) 

•RETORN 

•END 



Runstream TR DXDV 5 


This runstream implements the overall central difference method using either 
fixed or updated basis vectors. 

* 

t (TE DXDV 5) - CALCULATES DERIVATIVES OF TRANSIENT RESPONSE 
t USING TWO POINT CENTRAL DIFFERENCE OPERATOR 

t WITH UPDATED OR FIXED NODES 

* 

•XQT U1 

•REGISTER ST0RE(1 DXDV REGISTERS 11) 

•REGISTER RETRIEVE (1 DXDV REGISTERS 1 1) 

•RGI 

FDCH = .001 
FDNCH = .0001 
XLIB = S 
TLIB = 6 
OPT = 0 
RLIB = 14 
DRKETH0D=0 
DXMD=UPDATED 

•DC ALL, OPT (DXDV PARAMETERS) 

•SHOW 

I 

I INITIALIZATION FOR DERIVATIVE CALCULATIONS 

t 

•DCALL(TR.DPREP) 

$ 

$ LOOP OVER ALL DESIGN VARIABLES 

t 

! NDV = T0C,NJ(1 X ADS 1 1) 

! NCNT = NDV 
! IDV = 1 
•LABEL 10 
•iqT U1 

! IFLG = DS , 1 ,"IDV" ,1(1 XFLG ADS 1 1) 

*JZ(IFLG,100) 

•XQT U3 
RP2 

FORMAT 1 ’ ( 1H1 ,20X .41HBEGINNING SENSITIVITY CALCULATION. DV NO., 13) 

PRINT (1) "IDV" 

•LIBS "XLIB" 2 3 4 1 6 7 8 9 10 11 12 13 14 15 18 17 18 19 20 
•XQT AUS 
! IS = 1 
! NST = 2 
! SIGN =1.0 
$ 

I DO ANALYSIS FOR BOTH POSITIVE AND NEGATIVE STEPS 
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* 

•LABEL IS 
•XQT DCU 

COPT "XLIB" 1 XNAME ADS 1 1 
•XQT ADS 

DEFINE I = "XLIB" X ADS 1 1 
TABLE(NI=1,NJ="NDV") : X ADS 1 1 
TRAN(SOURCE=X) 

$ DIFFERENCE APPROPRIATE DESIGN VARIABLE 
! X = DS,1,"IDV",1("XLIB" X ADS 1 i) 

• DX = FDCH*X 

•IFC'DI" GT "FDCHM") : *GOTO 20 
! DX = FDNCH 
•LABEL 20 
! X = DX*SIGN + X 
TABLE, D : X ADS 1 1 
OPER = ISDN 
J="IDV" : "X" 

•C ALL ( SENS, DVDP) 

$ FORM PERTURBED MODEL 

*0NLINE=0 

•DCALL(HODEL) 

•0NLINE=1 
•xqT DCD 

COPT "XLIB" 1 TIME 
COPT "XLIB" 1 CA 
COPT "XLIB" 1 DMPD 
•IFC'DXMD" NE FIXE): *GOTO 30 
COPT "XLIB" 1 VIBR MODE 
*DCALL(TR,DIAG) 

•GO TO 40 
•LABEL 30 

•IFC’DXMD" NE UPDA): *GOTO 40 
•DCALL(TR, VECTORS) 

•LABEL 40 
•XqT AUS 

DEFINE X = VIBR MODE 1 1 1 , "NMDDES" 

DEFINE F * APPL FORC 1 
XTF = ITT(X,F) 

•DCALL(TR.REDM) VLIB=1 
•XqT DRX 

DTEI(DT="DT",METHOD="DRMETHOD",NTERMS="NTERMS") 

TR1(QILIB=1 , T2="T2" , qXlLIB=l , QI2LIB= 1 , LB="BLKSIZE" ) 
•DCALL(TR,DBACK 1) 

! SIGN = -1.0 

•LIBS "TLIB" 2 3 4 1 S 7 8 0 10 11 12 13 14 15 16 17 18 19 20 
*JGZ,-1(NST,15) 

* 

t COMPUTE DERIVATIVES USING CENTRAL DIFFERENCE OPERATOR 

I 

! TWDX = 2.0*DX 



! OVDX * l.O/TVDX 
! MOVD = - OVDX 
■1=1: ! M = NBCK 

•LIBS 1 2 3 4 S 6 7 6 9 10 11 12 13 14 IS 18 17 18 19 20 
•LABEL 80 
•XQT AOS 

! MM = DS, "I", 1,1(1 BACK LIST 1 1) 

! IEEE » TOC, IEEE (1 SEL "MM" MASK MASK) 

* JMZ ( IEEE , 90 ) 

DEFIME CPI = "XLIB" CEPT "MM" 

DEFIME CPO = "TLIB" CEPT "MM" 

DXDV "MM" "IDV" * SOMC'OVDX" CPI "MOVD" CPO) 

• XQT DCO 

PEIMT 1 DXDV "NM" "IDV" 

•LABEL 90 
! I ■ I ♦ 1 
•JGZ,-1(H,80) 

EEASE "XLIB" 

EEASE "TLIB" 

•LABEL 100 
! IDV = IDV + 1 
•JGZ,-l(MCNT, 10) 

$ 

•XQT 01 

•EEGISTEE EETEIEVEd DXDV EEGISTEES 1 1) 

•EETOEN 

•END 



Runstream TR DXDV 6 


This runstream implements the semi-analytical method with non-zero d$/dx. 
The called procedure TR DPHI determines how the basis vector derivatives are cal- 
culated. 

I 

$ (TR DXDV 6) - CALCULATES DERIVATIVES OF TRANSIENT RESPONSE 
I SEMI- ANALYTICALLY BUT WITH THE EFFECT OF CHANGING 

t NODES INCLUDED 

I 

•IQT U1 

* REGISTER ST0RE(1 DXDV REGISTERS 1 1) 

•REGISTER RETRIEVE(1 DXDV REGISTERS 1 1) 

•RGI 

FDCH = .001 
FDMCH = .0001 
XLIB » & 

OPT = 0 
RLIB = 14 
DRMETH0D=0 

*DCALL,OPT(DXDV PARAMETERS) 

•SHOW 

$ 

t INITIALIZATION FOR DERIVATIVE CALCULATIONS 

$ 

•DCALL(TR.DPREP) 

•XQT U1 

$ 

t LOOP OVER ALL DESIGN VARIABLES 

t 

! NDV * T0C,NJ(1 X ADS 1 1) 

! NCNT = NDV 
! IDV = 1 
•LABEL 10 

•LIBS "XLIB" 2 3 4 1 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
•XQT U1 

! IFLG = DS,1. "IDV", 1 ("XLIB" XFLG ADS 1 1) 

•JZClFLG , 100) 

•XQT DCU 

COPY "XLIB" 1 INANE ADS 1 1 
•XQT U3 
RP2 

FORMAT 1 ’ ( 1H1 , 20X , 41HBEGINNING SENSITIVITY CALCULATION. DV NO., 13) 

PRINT d) "IDV" 

•XQT AUS 

DEFINE X = "XLIB" X ADS 1 1 
TABLE(NI=1,NJ="NDV") : X ADS 1 1 
J=1,"NDV" : 1.0 


172 



J="IDV" : "FDCH" 

TRAN(SOURCE=X, OPERATION=MULT) 

$ CHECK FOK TOO SHALL A STEP 
! X = DS ,1 ,"IDV" , 1("XLIB" X ADS 1 1) 

! DX = FDCH*X 

*IF("DX" GT "FDCHM") : *GOTO 20 

• DX * FDMCH 

! I * X ♦ FDMCH 
TABLE.U : X ADS 1 1 
OPER = ISOM 
J="IDV" : "X" 

•LABEL 20 
*CALL(SENS,DVUP) 

! OVDX = 1.0/DX 
! HOYD = -OYDX 
$ FORK PERTURBED MODEL 
*0NLINE=0 
•DCALL(MODEL) 

$ CALCULATE DERITATITES OF MODES SHAPES 
•XQT AUS 

DEFINE KO = "XLIB" K SPAR 

DEFINE K1 = K SPAR 

DEFINE MO = "XLIB" "MNAME" 

DEFINE Ml = "MNAME" 

DKDY = SUM("OYDI" K1 "HOYD" KO) 

DMDV = SUMC'OYDI" Ml "MOYD" MO) 

•DCALL(TR,DPHI,3) 

•XQT DCU 

COPT "XLIB" 1 DT AUS 
COPT "XLIB" 1 DTEX AUS 

COPT "XLIB" 1 DCON AUS $ CONSTANTS FOR NEHMARK METHOD 
•XQT AUS 

DEFINE XO = "XLIB" YIBR MODE 1 1 1,"NH0DES" 

DEFINE KO * "XLIB" K SPAR 

DEFINE K1 = K SPAR 

DEFINE MO = "XLIB" "MNAME" 

DEFINE Ml = "MNAME" 

DEFINE DO = "XLIB" DAMP SPAR 

DEFINE D1 = DAMP SPAR 

DEFINE FO = "XLIB" APPL FORC 1 1 

DEFINE FI = APPL FORC 1 1 

DEFINE DXDY = DIDY AUS "IDY" 

$ 

t CALCULATE DERIVATIVE TERMS INVOLVING THE STIFFNESS MATRIX 

* 

DKX1 = PRQD(DKDV , XO) 

DKX2 = PROD(KO.XO) 

DKX3 = PROD (KO , DXDY) 

XDK1 = XTT(XO.DKIl) 

XDK2 = XTT(DXDV,DKX2) 

XDK3 = XTT(X0,DKX3) 
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XTMP = SUM ( XDK 1 , XDK2 ) 

XDKX = SUM ( XTMP, XDK3) 

$ 

$ CALCULATE DERIVATIVE TERMS INVOLVING THE MASS MATRIX 

t 

DMX1 = PROD(DMDV.XO) 

DMX2 = PROD(MO.XO) 

DMX3 = PROD (MO, DID V) 

XDM1 = XTY(XO.DMXi) 

XDM2 = XTT ( DXD V , DMX2 ) 

XDM3 = XTT(X0,DMX3) 

XTMP = SUM ( XDM1 , XDM2 ) 

XDMX = SQM ( XTMP, XDM3) 

* 

$ CALCULATE DERIVATIVE TERMS INVOLVING THE DAMPING MATRIX 

t 

! IDSP = TOC , IERR( 1 DAMP SPAR MASK MASK) 

•IFC’IDSP" NE 0): *GOTO 30 
DDDV » S0M("OVDX" D1 "MOVD" DO) 

DDX1 » PROD (DDDV, 10) 

DDX2 - PROD(DO.XO) 

DDX3 = PROD (DO, DXD V) 

XDD1 * ITT(I0,DDX1) 

XDD2 = XTY(DXDV , DDX2) 

XDD3 s XTT(X0,DDX3) 

ITMP = SUM ( XDD 1 , XDD2 ) 

IDDI = SUM ( XTMP, XDD3) 

•LABEL 30 

t 

$ CALCULATE DERIVATIVE TERMS INVOLVING THE FORCE VECTOR 

$ 

DFDV = SUMC’OVDX" FI "MOVD" FO) 

XF1 = XTT ( DXD V, FO) 

XF2 = XTY(XO , DFDV) 

XTF AUS = SUM(IF1, XF2) 

•XQT DCU 
PRINT 1 XTF AUS 
COPT "XLIB" 1 TIME AUS 
COPT "XLIB" 1 CA AUS 

$ 

t FORM THE RIGHT-HAND-SIDE PSEUDO LOAD VECTOR 

t 

•XQT DRX 

BACK(LRZ="BLKSIZE" ,PRINT=0) 

T = -1.0 XDMX s T = "XLIB" qX2 AUS 
•IFC’IDSP" EQ 0): T = -1.0 XDDX : T = "XLIB" QX1 AUS 
T = -1.0 XDKX : T = "XLIB" qi AUS 
Z = FH AUS 
•XQT DRX 

TRl(qXLIB=l,qXlLIB=l,qX2LIB=l,T2="T2",FHLIB=l,LB="BLKSIZE") 
•XQT DCU 



TOC 1 

* 

I BACK TRANSFORM FOR NECESSARY SENSITIVITIES OF PHYSICAL 
I QUANTITIES 

♦ 

♦DCALL(TR,DBACK,3) 

*XQTC DCU 
ERASE 1 
♦LABEL 100 
! IDV = ID? ♦ 1 
♦JGZ,-1(NCNT,10) 

I 

♦LIBS 1 2 3 4 6 6 7 8 9 10 11 12 13 14 15 18 IT 18 10 20 

*0NLINE=1 

♦xqT U1 

♦REGISTER RETRIE?E(1 DID? REGISTERS 1 1) 

♦RETURN 

♦END 
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Runstream TR DPREP 


This is a utility runstream used by all the sensitivity calculation runstreams. 
Its main task is to locate the critical points for all required response quantities. 

* 

* (TR, DPREP) - PREPARATION FOR SENSITIVITY CALCULATIONS 

t 

I 

t FORM CRITICAL POINT TABLES FOR RESPONSE QUANTITIES 

* 

! NBCK = TOC.NlC'qLIB" BACK LIST 1 1) 

! I = 1 : ! N = NBCK 
*LABEL 10 

! NM = DS , ”1" , 1 , 1 ("QLIB" BACK LIST 11) 

! I ERR = TOC, IERR( "QLIB" SEL "NM" MASK MASK) 

*JNZ(IERR,20) 

*XQT U10 

CRIT(Y="QLIB" HIST "NM" ,DT="DT" , NCRIT="NCRIT" , A 
CRPT="QLIB» CRPT "NM",CRTI="QLIB" CRTI "NM", k 
PCH= . 25) 

*XQT DCU 

PRINT "QLIB" CRTI "NM" 

PRINT "QLIB" CRPT "NM" 

•LABEL 20 
! I = I ♦ 1 
*JGZ,-1(N, 10) 

•XQT U1 

*(E4 PARAMETERS) EOFX 

♦RESET NFCT="NMODES" , NLIM="NMODES" 

RESET NIF="NMODES" 

IFS0URCE= "ILIB" VIBR MODE 1 1 
•EOFI 
•END 
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Runstream TR DBACK 1 


The TR DBACK n runstreams implement the different procedures for recovering 
the physical sensitivities. They all rely heavily on runstream TR CRPT which recovers 
a specific physical quantity at the critical points. Runstream TR DBACK 1 recovers 
the sensitivities using the mode displacement method and is used in the overall 
finite difference procedures. 

$ 

$ (TR DBACK 1) - BACK TRANSFORMATION FOR DERIVATIVES USING MODE 
$ DISPLACEMENT METHOD 

* 

! IIB = 1 : ! NNB = NBCK 
•LABEL 10 
•XQT AUS 

DEFINE X = VIBR MODE 1 1 1,"NM0DES" 

! NM = DS ,"IIB" , 1 , i("XLIB" BACK LIST 1 1) 

! I ERR = TOC , IERR ( "ILIB " SEL "NM" MASK MASK) 

•JNZ(IERR,200) 

* 

t DISPLACEMENTS 

* 

*IF("NM" NE DISP): *G0T0 30 
DEFINE IDJK = "ILIB" SEL DISP 
THAT VMOD = SVTRAN(IDJK.X) 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=1 LIB3=1 TNAME=VM0D CNAME=DISP Q=qX 
•GOTO 200 
•LABEL 30 

$ 

$ VELOCITIES 

$ 

*IF("NH" NE VELO): *G0T0 50 
DEFINE IDJK = "XLIB" SEL VELO 
THAT VVEL = SVTRAN(IDJK.X) 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=1 LIB3=1 TNAME=VVEL CNAME=VELO q=QXl 
•GOTO 200 
•LABEL 50 

$ 

$ ACCELERATIONS 

I 

*IF("NM" NE ACCE): *G0T0 70 
DEFINE IDJK = "XLIB" SEL ACCE 
THAT VACC = SVTRAN ( IDJK , X ) 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=1 LIB3=1 TNAME=VACC CNAME=ACCE q=qX2 
•GOTO 200 
•LABEL 70 

$ 
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* REACTIONS 

$ 

*IF("NH" ME EE AC) : *GOTO 90 
•GOTO 200 
•LABEL 90 

t 

$ STRESSES 

* 

*IF("NM" NE STRE): *GOTO 110 
•XgT ES 
RESET OPER=T 
IDQ = " XLIB" SEL STRESS 
0 = VIBR MODE 1 1 1 , "NMODES" 

T » THAT TSTRE 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=1 LIB3=1 TNAME=VSTRE CNAME=STRE Q=qX 

•GOTO 200 

•LABEL 110 

•LABEL 200 

! IIB = IIB ♦ 1 

•JGZ,-1(NNB,10) 

! IIB = FREE() : ! NMB = FREE() : ! LIB1 = FREE() : ! LIB2 = FREE() 

! LIB 3 = FREE() : ! TNAME = FREE ( ) : ! CNAME = FREE() 

! g = FREE() 

•END 



Runstream TR DBACK 2 


This runstream recovers sensitivities in the fixed mode, mode displacement 
version of the semi-analytical method. 

t 

$ (TR DBACK 2) - BACK TRANSFORMATION FOR DERIVATIVES USING 
t SEMI- ANALYTICAL METHOD 

$ UPDATE HISTORT 

* 

I 6/22/88 WHG - MODIFIED FOR SEMI-ANALYTICAL METHOD 
$ 6/21/88 WHG - CREATED FROM (TR, DBACK, 1) FOR UPDATED MODES 

* 

! I IB = 1 : ! NNB = NBCK 
♦LABEL 10 
♦XQT AUS 

! NM = DS,"IIB" , 1 , lC'XLIB" BACK LIST 1 1) 

! I ERR = TOC , IERR("XLIB" SEL "NM" MASK MASK) 

* JNZ ( IERR , 300 ) 

* 

t DISPLACEMENTS 

• 

*IF("NM" NE DISP): *GOTO 30 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2*"XLIB" LIB3=1 TNAME=VMOD CNAME=DISP Q=QX 
♦GOTO 200 
♦LABEL 30 

I 

$ VELOCITIES 

$ 

♦IF("NM" NE VELO): ♦GOTO 50 

♦DCALL(TR.CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=VVEL CNAME=VELO Q=QX1 
♦GOTO 200 
•LABEL 50 

* 

t ACCELERATIONS 

I 

•IF("NM" NE ACCE): •GOTO TO 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=VACC CNAME=ACCE Q=QX2 
♦GOTO 200 
•LABEL 70 

$ 

t REACTIONS 

t 

*IF("NM" NE REAC): •GOTO 00 
•GOTO 200 
♦LABEL 90 

t 

t STRESSES 

$ 
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*IF("NH" NE STEE): *GOTO 110 
$ FORM [S] [DQ/DV] 

♦DCALL(TR.CRPT) LIB1="XLIB" LIB2=”XLIB» LIB3=1 TNAME=VSTRE CNAME=STRE q=QX 
*XqT DCU 

CHANGE 1 CEPT STEE 1 1 CEPT STE1 1 1 
t FORM [DS/DV] [q] 

♦XqT ES 
RESET OPER=T 
IDq ■ "XLIB" SEL STRESS 
0 = "XLIB" VIBE NODES 1 1 1,"NM0DES” 

T = THAT YSTRE 
*XqT AOS 

DEFINE SO = "XLIB" THAT VSTRESS 

DEFINE SI = THAT VSTRESS 

THAT DSDV * SOH("OVDX" SI "HOVD" SO) 

♦DC ALL (TE, CEPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNANE=DSDV CNAHE=STEE q=qX 
♦XqT AOS 

DEFINE STR1 = CEPT STR1 

DEFINE STE2 = CEPT STEE 

CEPT STEE = S0H(STR1,STR2) 

♦GOTO 200 
♦LABEL 110 
♦GOTO 300 
♦LABEL 200 
♦XqT DCO 

CHANGE 1 CEPT "NM" 1 1 DXDV "NM" "IDV" 1 
COPT 1 "XLIB" DXDV "NM" "IDV" 1 
PRINT "XLIB" DXDV "NN" "IDV" 1 
♦LABEL 300 
! IIB = IIB ♦ 1 
♦JGZ,-1(NNB,10) 

■ IIB = FEEE() : ! NNB = FREE() : ! LIB1 = FEEE() : ! LIB2 = FREE ( ) 

! LIB3 = FREE() : ! TNAME = FEEE() : ! CNAME = FREE ( ) 

! q = FREE ( ) 

♦END 



Runstream TR DBACK 3 


This runstream recovers sensitivities in the semi- analytical method with non- 
zero d&/dx. 

* 

$ (TR DBACK 3) - BACK TRANSFORMATION FOR DERIVATIVES USING 
I SEMI -ANALYTICAL METHOD VITH CHANGING MODES 

* UPDATE HISTORY 

* 

* 8/22/88 WHG - MODIFIED FOR SEMI -ANALYTICAL METHOD 

I 8/21/88 WHG - CREATED FROM (TR, DBACK , 1) FOR UPDATED MODES 

| 

! IIB = 1 : ! NNB = NBCK 
•LABEL 10 
•XqT AUS 

DEFINE DX = DXDV AUS "IDV" 1 1,”NM0DES" 

! NM = DS ,"IIB" , 1 , 1("XLIB" BACK LIST 1 1) 

! I ERR = TOC , IERR("XLIB" SEL "NM" MASK MASK) 

•JNZ(IERR,300) 

t 

t DISPLACEMENTS 

* 

•IFC’NM" NE DISP): *G0T0 30 
DEFINE IDJK = "XLIB" SEL DISP 
TMAT DVMX = SVTRAN(IDJK ,DX) 

•DCALL(TR.CRPT) LIB1="XLIB» LIB2="XLIB" LIB3=1 TNAME=VMOD CNAME=DISP Q=qX 
•XqT DCU 

CHANGE 1 CRPT DISP 1 1 CRPT DSP1 1 1 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DVMX CNAME=DISP Q=qX 
•XqT AUS 

DEFINE D1 = CRPT DSP1 
DEFINE D2 = CRPT DISP 
CRPT DISP = SUM(D1,D2) 

•GOTO 200 
•LABEL 30 
* 

$ VELOCITIES 

$ 

•IF("NM" NE VELO): *GOTO 60 
DEFINE IDJK * "XLIB" SEL VELO 
THAT DVMX = SVTRAN( IDJK ,DI) 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2="XLIB'* LIB3=1 TNAHE=VVEL CNAME=VELO q=qXl 
•XqT DCU 

CHANGE 1 CRPT VELO 1 1 CRPT VEL1 1 1 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAHE=DVMX CNAME=VELO q=qXl 
•XqT AUS 

DEFINE VI = CRPT VEL1 
DEFINE V2 = CRPT VELO 
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CH.PT YELO = SDM(¥1,Y2) 

•GOTO 200 
•LABEL 50 

* 

t ACCELERATIONS 

$ 

•IFC'NN" NE ACCE): *GOTO 70 
DEFINE IDJK = "ILIB" SEL ACCE 
TKAT DVMX = SVTRAN(IDJK,DI) 

•DC ALL (TR, CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=¥ACC CNAHE=ACCE Q=qi2 
•XQT DCD 

CHANGE 1 CRPT ACCE 1 1 C&PT ACC1 1 1 

•DCALL(TR.CRPT) LIB1="XLIB» LIB2=1 LIB3="ILIB" TNAKE=DVHX CNAME=ACCE Q=QX2 
•XqT ADS 

DEFINE A1 = CHPT ACC1 
DEFINE A2 = C&PT ACCE 
C&PT ACCE = SDM(A1,A2) 

•GOTO 200 
•LABEL 70 
$ 

$ REACTIONS 

* 

•IFC'NM" NE &EAC) : *GOTO 80 
•GOTO 200 
•LABEL 90 

$ 

$ STRESSES 

$ 

•IFC'NN" NE STRE): *GOTO 110 
I FORM [S] [Dq/D¥] 

•DC ALL (TR, CRPT) LIB1="ILIB" LIB2="XLIB" LIB3=1 TNAHE=¥STRE CNAME=STRE q=qX 
•XqT DCU 

CHANGE 1 CRPT STRE 1 1 C&PT ST&l 1 1 
* FORM [DS/D¥] [q] 

•XqT ES 
RESET OPER=T 
IDq = "XLIB" SEL STRESS 
D = "ILIB" YIBR MODES 1 1 1,"NM0DES" 

T = TMAT YSTRE 
•XqT ADS 

DEFINE SO = "XLIB" TMAT YSTRESS 
DEFINE SI = TMAT YSTRESS 

TMAT DSDY = SDM("0¥DX" SI "MOYD" SO) 

•DCALLCT&.CRPT) LIB1="XLIB" LIB2=1 LIB3="ILIB" TNAME=DSDV CNAME=STRE q=qi 
•XqT DCD 

CHANGE 1 CRPT STRE 1 1 CRPT STR2 1 1 
$ FORM S [D PHI / DY] [q] 

•XqT ES 
RESET OPER=T 
IDq = "ILIB" SEL STRESS 
D = DXDY ADS "IDY" 1 1,"NH0DES" 
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T = THAT DSTRE 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAHE= DSTRE CNAME=STRE Q=QX 
•XQT ADS 

DEFINE STR1 = CRPT STR1 
DEFINE STR2 = CRPT STR2 
DEFINE STR3 = CRPT STRE 
TUP = SUH(STR1 ,STR2) 

CRPT STRE = SDM(THP , STR3 ) 

•GOTO 200 
•LABEL 110 
•GOTO 300 
•LABEL 200 
•XQT DCD 

CHANGE 1 CRPT "NH" 1 1 DXDV "NM" "ID?" 1 
COPT 1 "XLIB" DID? "NM" "IDV" 1 
PRINT "XLIB" DID? "NM" "ID?" 1 
•LABEL 300 
! IIB = IIB + 1 
•JGZ,-1(NNB,10) 

! IIB = FREE() : ! NNB = FREE( ) : ! LIB1 = FREE() : ! LIB2 = FREE() 

! LIB3 = FREE() : ! TNAME = FREE() : ! CNAME = FREE() 

! Q = FREE() 

•END 



Runstream TR DBACK 4 


This runstream recovers sensitivities in the fixed mode semi- analytical method 
using the mode acceleration method. 


(TE DBACK 4) - BACK TRANSFORMATION FOR DERIVATIVES USING 
SEMI- ANALYTICAL METHOD 
WITH THE MODE ACCELERATION METHOD 


•XqT AUS 

DEFINE X = "XLIB" VIBR MODE 1 1 1,"NM0DES" 

DEFINE E = "XLIB" VIBR EVAL 1 1 
DEFINE DKX ■ DKX AUS 
DEFINE DMX = DMX AUS 
DEFINE ROMG = "XLIB" ROMG AUS 
DEFINE XTDX = "XLIB" XTDX 
DEFINE XOMD = "XLIB" XOMD 
DEFINE XOME = "XLIB" XOME 

*IF("IDSP" NE 0): TABLE(NI="NMODES" ,NJ="NMODES") : RDDX 
$ CALCULATE VELOCITY TERM 
X0N1 = CBR(XOME,RDKX) 

XOMK = CBD(IOMI.ROHG) 

! NJDM = TOC,NJ("XLIB" XTDX MASK MASK MASK) 

! NBDM = TOC, NINJ( "XLIB" XTDX MASK MASK MASK) 

*IF("NJDM" NE "NBDM"): XOKC = CBR(XOMK,XTDX) 

•IF ("NJDM" Eq "NBDM"): XOKC = CBD (XOMK , XTDX) 

XOMC = CBR(XOME,RDDX) 

XqD = SUM(XOKC, -1.0 XOMC) 

# CALCULATE ACCELERATION TERM 
DKXO = CBD (DKX, ROMG) 

APPL FORC 887 = SUM(DKXO, -1.0 DMX) 

8 CALCULATE DERIVATIVE OF THE PSEUDO-STATIC TERM 
DEFINE USTAT = "XLIB" STAT DISP 1 1 
FSL1 = PROD(DKDV, USTAT) 

APPL FORC 888 = SUM(DFDV ,-l .0 FSL1) 

•XqT SSOL 

RESET SET=887, KLIB="XLIB" , KILIB="XLIB" , REAC=0 
•XqT SSOL 

RESET SET=888, KLIB="XLIB", KILIB="XLIB", REAC=0 

8 

8 LOOP OVER ALL RESPONSE qUANTITY TYPES 

8 

! IIB = 1 : ! NNB = NBCK 
•LABEL 10 
•XqT AUS 

! NM = DS ,"IIB" , 1 , 1 ("XLIB" BACK LIST 1 1) 

! IERR = TOC, IERR( "XLIB" SEL "NM" MASK MASK) 
*JNZ(IERR,300) 
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8 

8 DISPLACEMENTS 

8 

*IF("NM" NE DISP): *GOTO 30 
DEFINE IDJK = "ILIB" SEL DISP 
DEFINE DiCl = STAT DISP 887 
DEFINE DUST = STAT DISP 888 
DEFINE XTDI = "ILIB" XTDX 
DEFINE DACC = "ILIB" THAT DACC 
TMAT DUST = SVTRAN(IDJK .DUST) 

THAT DAC1 = SVTRAN ( ID J K , D AC 1 ) 

8 TELOCITY TERMS 
TMAT DTL1 = SVTRAN(IDJK ,XQD) 

$ 

*DCALL(TR,CRPT) LIB1="XLIB" LIB2=1 LIB3=»XLIB" TNAME=DUST CNAME=DISP Q=A 
T0CC(1 CRPT DISP 11): N2=DSP1 

*DCALL(TR,CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DVL1 CNAME=DISP Q=QX1 
T0CC(1 CRPT DISP 1 1) : N2=DSP2 

•DC ALL (TR, CRPT) LIB1="ILIB» LIB2="XLIB" LIB3=1 TNAME=DVEL CNAME=DISP Q=QX1 
T0CC(1 CRPT DISP 11): N2=DSP3 

•DCALL(TA.CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DAC1 CNAME=DISP Q=QX2 
T0CC(1 CRPT DISP 1 1) : N2=DSP4 

•DC ALL (TR, CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=DACC CNAME=DISP Q=QX2 
T0CC(1 CRPT DISP 1 1) : N2=DSP5 
•XQT U4 
TU 

DEFINE D1 = CRPT DSP1 
DEFINE D2 = CRPT DSP2 
DEFINE D3 = CRPT DSP3 
DEFINE D4 = CRPT DSP4 
DEFINE D5 = CRPT DSPS 

CRPT DISP = SUM(D1, D2, -1.0*D3, D4, -1.0*D5) 

•GOTO 200 
•LABEL 30 

$ 

8 VELOCITIES 

$ 

•IF("NM" NE VELO): *0010 SO 

•DC ALL (TR, CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=VVEL CNAME=VELO Q=QX1 
•GOTO 200 
•LABEL SO 

t 

t ACCELERATIONS 

$ 

•IF("NM" NE ACCE): *GOTO 70 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=VACC CNAME=ACCE q=QX2 
•GOTO 200 
•LABEL 70 

8 

8 REACTIONS 

8 
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*IF("NM" HE RE AC) : *GOTO 90 
•GOTO 200 
•LABEL 90 

t 

* STRESSES 

t 

*IF("NM" NE STRE): *GOTO 110 

* 

$ FORM [S] [DO/DT] 

* 

•LIBS 1 2 3 4 S 6 7 8 9 10 11 12 13 14 16 16 17 18 19 20 
•XQT ES 
RESET OPERsT 
XDQ * SEL STRESS 

U="XLIB" STAT DISP 888 : T * "XLIB" TMAT DTM1 

U="XLIB" XQD AOS : T * "XLIB" TMAT DTM2 

0="XLIB" STAT DISP 887 : T = "XLIB" TMAT DTM4 

•LIBS "XLIB" 2 3 4 1 6 7 8 9 10 11 12 13 14 16 16 17 18 19 20 
•XQT AUS 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DTM1 CNAME=STRE Q=A 
T0CC(1 CRPT STRE 1 1) : N2=STR1 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DTM2 CNAME=STRE Q=QX1 
T0CC(1 CRPT STRE 11): N2=STR2 

•DCALL(TR.CRPT) LIB1="XLIB" LIB2=»XLIB" LIB3=1 TNAME=SD CNAME=STRE Q=QX1 
T0CC(1 CRPT STRE 11): N2=STR3 

•DC ALL (TR, CRPT) LIB1="XLIB» LIB2=1 LIB3»"XLIB" TNAHE=DTM4 CNAME=STRE Q=QX2 
T0CC(1 CRPT STRE 11): N2=STR4 

•DC ALL (TR, CRPT) LIB1="XLIB" LIB2="XLIB" LIB3=1 TNAME=SP CNAME=STRE Q=QX2 
T0CC(1 CRPT STRE 11): M2=STR6 

t 

t FORM [DS/DV] [0] 

* 

•XQT ES 
RESET OPER=T 
IDq = "XLIB" SEL STRESS 
0 = "XLIB" STAT DISP 11 : T = TMAT SF 

0 = "XLIB" XOMD AUS : T = TMAT SD 

U = "XLIB” XOME AOS : T = TMAT SP 

•XqT AOS 

DEFINE SO = "XLIB" TMAT SF : DEFINE SI = TMAT SF 
TMAT DSF = SUM("OVDX" SI "MOVD" SO) 

DEFINE SO = "XLIB" TMAT SD : DEFINE SI = TMAT SD 

TMAT DSD = SOMC'OVDX" SI "MOVD" SO) 

DEFINE SO • "XLIB" TMAT SP : DEFINE SI = TMAT SP 

TMAT DSP = SUM("OVDX" SI "MOVD" SO) 

•DC ALL (TR, CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DSF CNAME=STRE q=A 
T0CC(1 CRPT STRE) : N2=STR6 

*DCALL(TR,CRPT) LIB1="XLIB" LIB2=1 LIB3="XLIB" TNAME=DSD CNAME=STRE q=qil 
T0CC(1 CRPT STRE) : N2=STR7 

•DCALL(TR.CRPT) LIB1="ILIB" LIB2=1 LIB3="XLIB" TNAME=DSP CNAME=STRE q=qX2 
T0CC(1 CRPT STRE) : N2=STR8 
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•XQT D4 
VO 


DEFINE 

SI 

S 

CRPT 

STR1 

DEFINE 

S2 

s 

CRPT 

STR2 

DEFINE 

S3 

= 

CRPT 

STR3 

DEFINE 

S4 

= 

CRPT 

STR4 

DEFINE 

SB 

= 

CRPT 

STR5 

DEFINE 

S6 

= 

CRPT 

STR6 

DEFINE 

S7 

= 

CRPT 

STR7 

DEFINE 

S8 

= 

CRPT 

STR8 

CRPT STRE 

= 

SUM(S1, S2, -1.0*S3 

•GOTO 200 




•LABEL : 

110 




•GOTO 300 




•LABEL 200 




•XQT DCU 





CHANGE 1 CRPT "NM" 1 1 DID? "NM" "IDV" 1 
COPT 1 "XLIB" DXDV "NM" "IDV" 1 
PRINT "XLIB" DXDV "NM" "IDV" 1 
•LABEL 300 


S6, -1.0»S7, -1 . 0*S8) 


! IIB = IIB + 1 
•JGZ,-1(NNB,10) 

! IIB = FREEO : ! NNB = FREE() : ! LIB1 = FREE() : ! LIB2 = FREE() 
! LIB3 = FREEO : ! TNAME • FREE() : • CNAME = FREEO 
! Q * FREEO 
•END 
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Runstream TR DPHI 3 


This runstream implements the modified modal method for calculating eigen- 
vector derivatives and is called from sensitivity calculation runstream TR DXDV 6. 

* 

$ (TR,DPHI ,3) - CALCULATE EIGENVECTOR DERIVATIVES USING THE 
$ MODIFIED NODAL METHOD 

* 

*iqT AUS 

INLIB=10 : 0UTLIB=10 
DEFINE MO = "ILIB" "MNAME" 

DEFINE DK - 1 DXDV SPAR 
DEFINE DM = 1 DMDV SPAR 

DEFINE XO = "XLIB" VIBR MODE 1 1 1,"NM0DES" 

DEFINE AJK = AJK 

DEFINE EV = "XLIB" VIBR EVAL 1 1 
MX - PROD (DM, XO) 

AKK = XTTD(-.B XO,MI) 

T ABLE ( NI= "NMODES " , N J a"NMODES " ) : AJK 

TABLE(NI=1 ,NJ="NM0DES") : UNIT : J=l, "NMODES" : 1.0 

! J = 1 : ! NJ = NMODES 
•LABEL 10 

! EJ * DS,"J», 1,1 ("XLIB" VIBR EVAL 1 1) 

! MEJ * -EJ 

DEFINE XJ * "XLIB" VIBR MODE 1 1 "J»,”J» 

DKDM a SUM(DK,"MEJ" DM) 

MX • PROD(DKDM.XJ) 

DLAM a ITT(XJ.MX) 

! DLAM a DS, 1,1, 1(10 DLAM AUS 1 1) 

AF1 a PROD ("DLAM" MO, XJ) 

AF2 a SUM( AF1 -1.0 MI) 

*IF("J" EQ 1): 11 APPL FORC = UNION (AF2) 

*IF("J" NE 1): 11 APPL FORC = UNION, U(AF2) 

AA a XTT(XO.MX) 

DENI a SUM(-1.0 EV, "EJ" UNIT) 

DEN2 a PROD (EV, DENI) 

TABLE, U : DEN2 : Ia"J" : Jal : 1.0 
FACT a RECIP(DEN2) 

AAB a PRODC’EJ" FACT.AA) 

DEI : OPER=XSUM : DEST,U=AJK AUS 

SOURCEaAAB : JS=1 : JD="J" : EX1 

SOURCEaAKK : IS="J" : JS=1 : IDa"J" : JDa"J" : Eli 

! J a j + l 

•JGZ ,-l(NJ , 10) 

DXDV a CBR(X0 , AJK) 

•XqT SSOL 

RESET KLIBa"XLIB" , KILIB="XLIB" , qLIB=ll, REACaO, EP=0 
•XqT AUS 
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DEFINE D1 = 10 DID? ADS 
DEFINE D2 * 11 STAT DISP 
DID? ADS "IDF" = SDM(D1 ,D2) 
•DELETE 10 
•DELETE 11 
•RETURN 
•END 



Runstream TR CRPT 


This is a utility runstream which performs the transformation from modal to 
physical basis for a single response quantity at a set of critical times. It is heavily 
used by the TR DBACK runst reams. 

A 

A (TR.CRPT) - FORK CRITICAL POINT RESPONSE TABLE 

A 

*0NLINE=0 
*XQT AUS 

! NCRIT * TOC ,NI ("LIBl" CRTI "CNAKE" 1 1) A NUMBER OF TIKES 

! ND * T0C,NI("LIB2" TKAT "TNAME" 1 1) A NUMBER OF RESP. QUANTITIES 

! NQ = T0C,NJ("LIB2" TKAT "TNAME" 1 1) A NUMBER OF MODES 

! NJQ = T0C,NJ("LIB3" "Q" AUS MASK MASK) 

! ISTP = 0 

A LOOP OVER ALL RESPONSE QUANTITIES 
! I = 1 
! N = ND 

INLIB = 21 : OUTLIB = 21 
•LABEL 20 
DEI 

S0URCE="LIB2" THAT "TNAME" 

ID * 1 : IS="I" 

DEST=T0NE "TNAME" "I" 1 
Ell 

TABLE(NI="NQ" , NJ="NCRIT") : XBAR CRIT "I" 

A LOOP OVER NUMBER OF CRITICAL POINTS 
! II = 1 
! NN = NCRIT 
•LABEL 40 
DEI 

! TIME = DS , "II" , 1 , ”I”("LIB1" CRTI "CNAME" 1 1) 

ISTP = TIME/DT + .6 
! ISTP * ISTP ♦ 1 
! IB = ISTP/NJQ 
! 1ST = IB*NJQ 

*IF("IST" NE "ISTP"): ! IB = IB ♦ 1 
! J = IB - 1 * NJQ : ! J = ISTP - J 

SOURCE = "LIB3" "Q" AUS MASK MASK "IB" : JS="J" : DEST,U=XBAR CRIT "I" 

JD = "II" : Ell 
! II = II ♦ 1 
*JGZ,-1(NN,40) 

DEFINE T = TONE "TNAME" "I" 

DEFINE XBAR = XBAR CRIT "I" 

Cl AUS "I" = RPROD(T.XBAR) 

T0CC(CI AUS "I") : NJ=1 
DEFINE Cl = Cl AUS "I" 

*IF("I" EQ 1): 1 CRPT "CNAME" = UNION(CI) 
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*IF("I" GT 1): 1 CRPT "CNAME" = UNION, U(CI) 

! I = I ♦ 1 
*JGZ,-1(N,20) 

*ONLINE=l 

! NC&IT = FREE() : ! ND = FEEE( ) : ! NQ = FREE ( ) ! NJq = FREE() 
! ISTP = FREE() : ! IB = FREEQ : ! 1ST = FREE() 

•END 


Runstream TR PI AG 


This is a utility runstream that solves the reduced order eigenproblem based 
on a given set of basis vectors to uncouple a reduced system. 

$ 

$(TR DIAG) - MAKE AN ARBITRARY VECTOR SET K AND K ORTHONORMAL 

t 

• XQT AUS 

0UTLIB=10: INLIB=10 

DEFINE K=1 K: DEFI M*1 "MNAME" 

DEFINE X = "VLIB" VIBR MODE 
IJC0DE=10000 

KX=PROD(K,I): STN X 10000 "NMODE" = XTTS(X.KX) 

MX=PR0D(M,X) : STN H 10000 "NMODE" = XTTS(I,MX) 

! ZER0=NM0DE- 1 

• JZ (ZERO, 1003) 

• XQT STRP 

RESET S0URCE=10, DEST=10 

• JGZ (ZERO, 1004) 

• LABEL 1003 

• XQT AUS 

0UTLIB=10: INLIB=10 
! K=DS 2 1 1(10 STN K MASK MASK) 

!M=DS 2 1 1(10 STN M MASK MASK) 

•EYALsK/M 

TABLE(NI=1 ,NJ=l) : STS EVEC: J=l: 1.0 
TABLE(NI=1,NJ=1): STS EVAL: J=l: "EVAL" 

• LABEL 1004 

• XQT AUS 

0UTLIB=10: INLIB=10 
DEFINE E^STS EVEC 
DEFI X = "VLIB" VIBR MODE 
X ORTH 1 1=CBR(X,E) 

DEFINE X=X ORTH 1 1 

"VLIB" VIBR MODE 1 1=UNI0N(X) 

DEFINE E=SYS EVAL 

"VLIB" VIBR EVAL 1 1=UNI0N(E) 

•END 
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